Meanphotonnumber dependent variational method to the Rabi model
Abstract
We present a meanphotonnumber dependent variational method, which works well in whole coupling regime if the photon energy is dominant over the spinflipping, 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 rotatingwave 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 meanphotonnumber dependent variational parameter. For the excited states,our method yields considerable improvements over the generalized rotatingwave approximation. The variational method proposed could be readily applied to the more complex models, for which an analytic formula is difficult to be formulated.
I Introduction
The Rabi model describes a twolevel system interacting with a singlemode 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 solidstate semiconductorGunter2004 (); Deveaud2007 (); Cristofolini2012 (); Carusotto2013 (); wen (), which provides a perfect test bed to explore the physics of lightmatter 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 multimodel Rabi model, i.e., the spinboson model Leggett_RMP (); weiss_book (), for which an exact solution is quite difficult to obtain, and a wellestablished 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 rotatingwave approximation (RWA) jc (). Working in the nearresonance and weakcoupling regime, the RWA neglects the counterrotating terms in the interaction and results in the JaynesCummings (JC) model jc (). It has served as a basic starting point in understanding many quantum phenomena involved in lightmatter 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 counterrotating term becomes important due to the strong circuit1 (); circuit2 () or the ultrastrong coupling ultra () between the bosonic field and the twolevel 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 twolevel system. After working in the displaced oscillator basis, it takes the frequency of the twolevel system as perturbation and results in a truncated Hamiltonian with the interaction effects collected in a renormalization factor to the frequency of the twolevel system. In 2007 grwa (), the AA was improved by considering the RWAtype 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 twolevel 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
(1) 
where and are the annihilation and creation operators of the quantized singlemode bosonic field with frequency , is the Pauli matrix for the twolevel system with level splitting , and are the transition operators between the two levels, and is the coupling strength. Here, for convenience of comparison we follow the notations in Ref.grwa () to use spinflipping for the levelsplitting term instead of commonly used in quantum optics scully (). However, these two notations can be transformed into each other by a rotation on the twolevel system. According to the tuning relationship between the twolevel system and the field, the model takes three cases: resonance (), positive detuning () and negative detuning (). Throughout the paper we take as unit of energy.
Essentially, the existing approximate methods can be formulated in two ways: One is to truncate Eq. (1) into JClike exactly solvable form, and the other is to expand Eq. (1) on a proper basis and then truncate the obtained matrix into the blockdiagonal 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 counterrotating terms in Eq. (1) yields the RWA Hamiltonian
(2) This is the JC Hamiltonian jc (), which is exactly solvable. Its eigen solution reads
(3) with the ground eigenenergy , which is just the JC energy ladder Fink2008 ().

AA: Performing a unitary transformation with to Eq. (1), one obtains with adiabatic ()
(4) Here , , and
with being the associated Laguerre polynomial (see Appendix A). In the small case[], keeping only the zeroth order term of and in is a good approximation, which leads to
(5) whose eigensolution can be evaluated readily as
(6) with being the eigenstates of and being the Fock state. After the inverse transformation, through representing the by the original basis, one gets the eigenstate under the AA:
(7) 
GRWA: Going beyond the AA, one further considers the zeroth order term in G() involving oneexcitation terms. Only considering the “energyconserving” oneexcitation terms, one arrives at the GRWA Hamiltonian grwa ()
(8) On the basis of , Eq. (8) is blockdiagonalized with subblocks
(9) which gives a pair of eigenvectors . The offdiagonal entries are defined by
(10) Thus, the eigenstates of Eq. (8) read as
(11) The states to the original Hamiltonian (1) are obtained by the inverse transformation:
(12) while the ground state is the same as that of AA.

GVM: Different from the above two methods, the parameter here is not fixed but is optimized by minimizing the groundstate energy GVMground ()
(13) which results in the equation to determine as . Since it cannot be solved analytically, Zhang et al. GVMground () took the following approximate solution
(14)
Below we address the conditions under which the above methods work well. The RWA is valid in the very weak coupling regime () and under the nearresonance () 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 or large , the term of displaced oscillator is dominant in (4) and the terms can be treated as perturbation. Thus the validity of the AA lies in strong coupling limit () or negative detuning () regime. Because the GRWA further keeps all oneexcitation “energyconserving” 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 groundstate 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 in the weak coupling regime, especially for the ground state. Third, it predicts an incorrect independent mean photon number due to the fixed . The GVM improves the accuracy of the groundstate energy and its mean photon number behavior captures the 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 offdiagonal 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 with , Eq. (1) can be rewritten as
(15) 
with and . Discarding the remote offdiagonal elements leads to a blockdiagonal matrix
(16) 
By diagonalizing Eq. (16), one finds the eigen solution
(17)  
(18) 
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
(19) 
with . Then dropping the remote offdiagonal matrix elements gives rise to
(20) 
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 photonnumber 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 blockdiagonalized 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
(21) 
where are the coherent states. For the AA/GRWA, , while for the GVM, is approximately given by Eq. (14). This motivates us to take Eq. (21) as our trial state but completely free the parameter . 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 spinflipping in the model. It is noticed that if the spinflipping term is neglected, the model reduces to a displaced oscillator Hamiltonian with two degenerated ground states and . Considering an infinitesimal spinflipping, the degeneracy is lifted and their linear combination, namely, is just the trial state Eq.(21) with . When increasing the spinflipping, the competition between the displaced oscillator and the spinflipping motivates us to free the as a variational parameter. In addition, is also eigenstate of parity operator , 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:
(22)  
(23) 
Obviously, the parameter is optimal if the projection is exactly equal to one, where is the exact ground state. Unfortunately, a simple expression of is unknown (though a series expression of can be given but it is quite useless due to the infinite series form). We here adopt an approximate but accurate enough form for . Note that a unitary transformation can recast to , which can be regarded as the zeroth order approximation for . Taking the complete orthogonal basis of into consideration, we can construct the perturbative corrections of . For perturbative calculation, we choose the basis of AA to be the complete orthogonal basis of . 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 to the first order of the perturbation as
(24) 
where , and with . Note that the primed summation excludes the ground state itself with the label , and will vanish if state has a different parity from . Then we can calculate
(25) 
For given and , if is minimized by choosing optimal , then the obtained 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 is much larger than , the twolevel splitting in Eq. (1) can be safely neglected. Then one can verify that Eq. (25) takes the form , which has an optimal value only when . 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 , which has the optimal value when . This recovers the result under the GVM in Ref. GVMground ().
For other cases, the analytic evaluation of the optimization on or is difficult. We should resort to numerical evaluation of the expression of , which has much less numerical work than that of exact numeric. Figure 1 shows as the function of in different and . The leading four components in the summation of are also plotted. We can see the following characters: (i) When is sufficiently large [see Fig. 1(a)], all the components except are negligible. Thus, is a good substitution of for the minimization. (ii) When is comparable to , for becomes sizable [see Fig. 1(b)]. However, still has only one minimum, which means that still can act as a substitution of for the minimization. (iii) When is small [see Fig. 1(c)], for become important and shows two minimums. Therefore, none of can be taken as a substitution of for the minimization. The two minimums of should be considered equally. (iv) For sufficiently small [see Fig. 1(d)], the series of lose convergency and a multiminimum structure of appears. This complicated structure indicates that the coherent state form of the trial wavefunction Eq. (21) cannot capture the physics dominated by the spinflipping and our scheme is no longer valid in this regime.
Figure 2 shows the optimal and the corresponding for the negativedetuning (), the resonance (), and the positivedetuning () cases. For the twominimum situation in the positivedetuning [see Fig. 2(a) and Fig. 2(d)] case, the optimal can be evaluated effectively as
(26) 
where and are the two minimum positions of , and and are their corresponding . We find that is dependent of the coupling strength, which is quite different from the fixed result under the AA/GRWA grwa (). Furthermore, approaches to in the weak coupling limit, which is consistent with the analytic result of obtained under the GVM. And approaches to in the strong coupling limit, which is consistent with the analytic result of obtained under the AA/GRWA. In the whole parameter range, shows little deviation from , which indicates that our obtained is almost the same as the exact ground state. For the oneminimum situation in the resonance [see Fig. 2(b) and (e)] and in the negativedetuning [see Fig. 2(c) and (f)] cases, where the optimal is determined by minimizing , it is interesting to find that the change scope of converges closer and closer to one. It means that our scheme performs better and better with the increase of the .
Figure 3 shows the groundstate energy and the mean photon number as a function of the coupling strength for different detuning cases evaluated by different methods. We stress that although 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 matches well with the exact one even in the whole coupling regime. The improvement of our scheme to becomes more obvious with the decrease of . 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 , the AA/GRWA and the GVM become worse and worse, but our results remain its good performance in evaluating and , 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 , the shows a multipleminimum 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 spinflipping becomes dominant.
iii.2 Applications to excited states
Our variational method can be also applied to the excited states. The GRWAform excited state is adopted to be the trial state, but with unfixed parameter
(27) 
This trial state possesses a definite parity. As in the ground state, the perturbation scheme is again employed to determine the optimal value of . For convenience, we still discuss in the transformed representation. The zeroth order Hamiltonian is , and the perturbation is . is determined by maximizing the projection , where is the exact excited state corresponding to . can be evaluated perturbatively as
(28) 
where and with . Here, similarly to (24), the primed summation excludes the trial state itself with the label . Then can be calculated by optimizing . The corresponding energy and mean photon number are
(29)  
(30)  
In the large limit, approaches , where
(31) 
Then the optimal takes , which recovers the result of the GRWA. In the small limit, approaches . Then the optimal reads . In both of the two limits, the chosen is independent of excitation label . It means the series of the approximate states hold the orthogonality. For other coupling cases, where can be extracted by simple numerics, one can expect that depends on excitation label . Exactly speaking, differences in 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. , 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 , 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 experimentrelated parameters in Ref.ultra (), which reads , and . This is a large detuning case and the detuning is also much larger than the coupling strength since and . 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 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 positiondisplaced 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 widelyused 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 multimode Rabi model, i.e., the socalled spinboson 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 transforms the model Hamiltonian into in the new representation. With the formula
(32) 
where if and otherwise, one can obtain
(33) 
The terms of and can be expanded in powers of and according to formula
(34) 
Below we will use the associated Laguerre function defined by
(35) 
and and the Laguerre function
(36) 
For the factor , one has
(37) 
where
(38) 
Here is particle number operator. Set .
(39)  
For ,
(40)  
For ,
(41)  
By the definition of the function
(42) 
one can expand Eq.(43) as
(43) 
By the same way, one has
(44) 
Substitute Eq.(43) and Eq.(44) into the transformed Hamiltonian Eq.(33), the expansion Eq.(4) is obtained.
References
 Rabi I I 1936 Phys. Rev. 49 324 Rabi I I 1937 Phys. Rev. 51 652
 Walther H, Varcoe B T H, Englert B and Becker T 2006 Rep. Prog. Phys. 69 1325
 Raimond J M, Brune M and Haroche S 2001 Rev. Mod. Phys. 73 565
 Holstein T. 1959 Ann. Phys. (Amsterdam, Neth.) 8 325
 Braak D 2011 Phys. Rev. Lett. 107 100401
 Solano E 2011 Physics 4 68
 Romero G, Ballester D, Wang Y M, Scarani V and Solano E 2012 Phys. Rev. Lett. 108 120501
 Restrepo J, Ciuti C, and Favero I 2014 Phys. Rev. Lett. 112 013601
 Crespi A, Longhi S, and Osellame R, 2012 Phys. Rev. Lett. 108 163601
 Wallraff A, Schuster D I, Blais A, Frunzio L, Huang R S, Majer J, Kumar S, Girvin S. M and Schoelkopf R. J 2004 Nature (London) 431 162
 Wallraff A, Schuster D I, Blais A, Frunzio L, Majer J, Devoret M H, Girvin S M and Schoelkopf R J 2005 Phys. Rev. Lett. 95 060501
 Niemczyk T, Deppe F, Huebl H, Menzel E P, Hocke F, Schwarz M J, GarciaRipoll J J, Zueco D, Hümmer T, Solano E, Marx A and Gross R 2010 Nature Phys. 6 772
 FornDiaz P, Lisenfeld J, Marcos D, GarciaRipoll J J, Solano E, Harmans C J P M, and Mooij J E 2010 Phy. Rev. Lett. 105 237001
 Günter G, Anappara A A, Hees J, Sell A, Biasiol G, Sorba L, De Liberato S, Ciuti C, Tredicucci A, Leitenstorfer A and Huber R 2009 Nature 458 178
 Deveaud B Ed. The Physics of Semiconductor Microcavities (WileyVCH, Weinheim, 2007).
 Cristofolini P, Christmann G, Tsintzos S, Deligeorgis G, Konstantinidis G, Hatzopoulos Z, Savvidis P, and Baumberg J 2012 Science 336 704
 Wen P, Christmann G, Baumberg J J and Nelson K A 2013 New J. Phys. 15 025005
 Carusotto I and Ciuti C 2013 Rev. Mod. Phys. 85 299
 Larson J 2012 Phys. Rev. Lett. 108 033601
 Caldeira A O and Leggett A J 1981 Phys. Rev. Lett. 46 211
 Weiss U Quantum Dissipative Systems, 3rd ed. (World Scientfic, Singapore, 2008)
 Casanova J, Romero G, Lizuain I, GarciaRipoll J J, and Solano E 2010 Phy. Rev. Lett. 105 263603
 Hausinger J and Grifoni M 2010 Phys. Rev. A 82 062320
 He S, Wang C, Chen Q H, Ren X Z, Liu T and Wang K L 2012 Phys. Rev. A 86 033837
 Chen Q H, Zhang Y Y, Liu T and Wang K L 2008 Phys. Rev. A 78 051801
 Jaynes E T and Cummings F W 1963 Proc. IEEE 51 89
 Shore B W and Knight P L 1993 J. Mod. Opt. 40 1195
 Mabuchi H and Doherty A C 2002 Science 298 1372
 Irish E K, GeaBanacloche J, Martin I, and Schwab K C 2005 Phys. Rev. B 72 195410
 Irish E K 2007 Phys. Rev. Lett. 99 173601
 Zhang Y, Chen G, Yu L, Liang Q, Liang J Q and Jia S 2011 Phys. Rev. A 83 065802
 Yu L, Zhu S, Liang Q, Chen G and Jia S 2012 Phys. Rev. A 86 015803
 Scully M O and Zubairy M S Quantum Optics, (Cambridge University Press, 1997).
 Fink J M, Göppl M, Baur M, Bianchetti R, Leek P J, Blais A, and Wallraff A 2008 Nature (London) 454 315
 Ashhab. S and Nori. F 2010 Phys. Rev. A 81 042311
 Vojta M 2012 Phys. Rev. B 85 115113
 Tong Q J, An J H, Luo H G and Oh C H 2011 Phys. Rev. B 84 174301