The quantum Rabi model^[1,2] plays a fundamental role in many fields of modern physics, such as quantum optics,^[3] quantum information,^[4] and condensed matter physics.^[5] It has been used to characterize a system composed of a two-level atom and a single mode bosonic field with a linear coupling in between, namely, the so-called cavity quantum electromagnetic dynamics (QED) system.^[3] Recently, the model has also been realized in many other promising quantum systems, e.g., the superconducting circuits system,^[6,7] trapped ion,^[8] and so on. A remarkable experimental achievement involves the realization of the so-called ultra-strong coupling.^[9–11] The realization from strong coupling to ultra-strong coupling and to even the deep-strong coupling inspires theorist to explore the interesting physics involved in the model. In fact, the quantum Rabi model has been extensively studied theoretically and a series of important features have been found, e.g. its integrability,^[12] geometric phase,^[13,14] and the polaron physics.^[15–17] In particular, more recently, quantum phase transitions and its universality class have been found in such a finite degree of freedom system.^[8,18–20]

The starting point of the theoretical research on such a quantum model lies in solving its eigenvalue problem. Although the form of the model looks quite simple, it is not easy to solve exactly. In fact, it is until 2011 that the model has been exactly solved by Braak.^[12] This pioneering work triggered a series of exact solutions in the related models.^[21–25] However, since the key technique of the exact solution involves solving the zero points of a transcendental function, using it to analyze the physical properties of the model is not trivial. Hence, people are still trying to find some intuitive analytic approximation to calculate and understand the model. The earliest and maybe the most famous approximation is the so-called rotating wave approximation (RWA). However, this approximation is only applicable to weak coupling with near-resonance condition, and thus it fails to satisfy the requirement of the experiment in current ultra-strong coupling regime.^[26] In 2007, Irish proposed the so-called general rotating wave approximation (GRWA),^[27] which does succeed in presenting the correct structure of the energy spectrum over a large parameter regime. It also provides a framework to analytically treat many models which are cousins of the quantum Rabi model.^[28–42]

We notice that the so-called bias term is usually involved in the experimental circuit QED system.^[26] Therefore, the study of the Rabi model with the bias term is indispensable. However, most of theoretical studies on the quantum Rabi model focus on the unbiased case. Although there already exists an analytical discussion based on the GRWA for the biased quantum Rabi model,^[30] we find that there is still space to improve the accuracy in the ultra-strong coupling regime, especially when the resonator frequency is small. Moreover, the ground state mean photon number and correlation predicted by GRWA are independent of the two-level energy splitting and the bias strength. It is actually not true. For these reasons, we show in this work an improvement to the GRWA by introducing an analytical variational method. Our method can be obtained by two equivalent ways, namely, the unitary transformation method and the trial ground state wavefunction method. The key point is that we determine the optimal value (by minimizing the ground state energy) of a variational parameter λ introduced in our method. To show the validity of our method and its improvement over the GRWA, we compare our results to the numerical exact ones and the ones from GRWA respectively. As an approximation, we also discuss its applicable condition.

The Hamiltonian of the biased quantum Rabi model reads^[43–46]

Our main idea of the variational method is to introduce an unfixed parameter λ in a unitary transformation (see below). The optimal value of λ can be determined by minimizing the ground state energy. There are two ways to develop the variational scheme. One is to perform a λ-dependent unitary transformation, the other is to directly introduce a λ-dependent ground state wave function. The two ways are proved to be equivalent to each other.

3.1. Variational method approached by a unitary transformation

We start our variational scheme by performing a unitary transformation similar to that in the GRWA.^[27,30] Thus the Hamiltonian in Eq. (2) is transformed into , where . Here λ is introduced as a variational parameter as in Ref. [34], which is in contrast to the GRWA where it is fixed to . We will see later that this difference allows us to improve the GRWA. The transformed Hamiltonian is explicitly written as

where The hyperbolic sine and cosine terms can be further expanded as^[47] The function F_m is defined aswhere m and n are integers, andis the Laguerre polynomial. Up to now, is exact, but it is still unsolvable. Therefore, we have to take some approximations to treat it analytically.

In this work, since we are interested in the low energy properties of the model, the photon excitations can be neglected in . Thus we approach an effective approximate Hamiltonian

The Hamiltonian in Eq. (11) is analytically solvable. We can choose the basis as , where are the eigenstates of with eigenvalues ±1, and is the Fock state. The effective Hamiltonian can be readily divided into invariant subspace as

where , and .

The eigenvalues of the matrix in Eq. (12) can be obtained as

and the corresponding eigenvectors arewhere

With the λ-dependent energy spectra and wave functions in hand, we can explore particularly the ground state properties and determine the optimal value of λ. Specifically, the ground state energy is

The λ can be determined by minimizing the ground state energy , namely, solving the equationIt can be solved in a numerical way or approximately aswhere .

In the limit cases, we can further simplify the formula of λ. In the large g limit, λ approaches to , which recovers the GRWA.^[27,30] In the small g case, there is

which is the same as λ₀. The small g case involves the regime where g is small but comparable to ω. This limit covers the typical ultra-strong coupling regime, which recently attracts many experimental interests and is our main focus in this work.

Having determined the value of λ, we fully obtain the analytical formula of the ground state wavefunction, and the physical observables can be thus calculated.

3.2. Variational method approached by a trial ground state wavefunction

Before going on, we present an alternative way to derive the variational scheme. Instead of performing the unitary transformation, we can in fact directly assume the form of trial ground state wave function. Remind that the trial wave function in the unbiased Rabi model takes the form

where is a coherent state. The equality of weight for up and down spins is due to the parity symmetry. Naturally, if the parity symmetry breaks down, the weights of up and down spins are no longer equal. Therefore, the trail ground state wave function can be assumed asFor normalization, we set and . Thus the ground state energy should beIn the last line, the equal condition is taken if and only if and . Thus the ground state energy iswhich recovers Eq. (17).

With the analytical expressed ground state wave function (Eq. (14) or Eq. (22)) in hand, we can readily calculate the physical observables. Since we are interested in the ultra-strong coupling regime, the λ is given by the analytical formula in Eq. (20). The quantum average of a physical observable Q is In this work, we take ground state energy , mean photon number , spin orientation , and correlation as examples. Besides the ground state energy in Eq. (17), the rest ones are

We first consider the ground state energy. In Fig. 1, we see that our variational result agrees with the numerical exact one very well. Comparing to the result from GRWA, we find that our variational method gives a lower energy curve. It means that in spirit of variational way, the method successfully achieves the energy minimum. From Fig. 1, we can also obtain some detail properties of the variational method. In Figs. 1(a) and 1(b), we see that the deviation between the GRWA and numerical one becomes large when the coupling strength increases. It can be understood that λ is proportional to g, and thus the increasing of g enlarges the difference between the optimal λ in our method and in the GRWA. For the same reason, from Fig. 1(c), we infer that λ is highly relevant to Ω, since the difference between our method and the GRWA becomes evident when Ω is large. Figure 1(d) shows the ground state energy as a function of the bias, from which one notes that our method gives an excellent result in comparison to the numerical exact one. In comparison to the GRWA, our improvement is definite.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) The ground state energy E0 for different parameters: (a) , Ω = 0.5, (b) , Ω = 5, (c) , g = 0.5, (d) g = 0.2, Ω = 5.

Since the mean photon and the correlation are simply depend only on λ, we can see the importance of the choice of λ in detail. Quantitatively, we see that our calculation agrees with the numerical exact result very well, while the GRWA presents an incorrect result. Qualitatively, our method seizes the important physics that the mean photon number and correlation are both dependent on Ω and (see Figs. 2(c)–2(f)). This feature is ignored by the GRWA. Our method is efficient in calculating most of the physics quantities, such as the average value of spin orientation (see Fig. 3).

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) The mean photon number and correlation for different parameters: (a), (b) , Ω = 5, (c), (d) , g = 0.1, (e), (f) g = 0.2, Ω = 0.1.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) The spin orientation σx for different parameters: (a) , Ω = 0.5, (b) , Ω = 5, (c) , g = 0.2, (d) g = 0.5, Ω = 0.1.

It would be valuable to briefly discuss the applicable regime of our method. The main idea of the GRWA is to consider the Ω term in the Hamiltonian as a perturbation. Thus, the GRWA would work well when Ω/ω is not large (see Fig. 4(a)). If Ω is large, the performance of the GRWA will depend on the value of g. Here we discuss three cases. Case I: g is large enough, the oscillator term and the coupling term together construct a displaced oscillator and they dominate. Hence the Ω term can still be considered as a perturbation and the GRWA works well in this case. Case II: g is small enough, the Ω term is non-perturbative, and the displacement parameter λ should be optimized. The λ optimized GRWA is our variational method. Case III: g is neither large nor small, we find that no matter how to optimize the λ, we can not lift the deviation between our result and the numerical exact one (see Fig. 4(b)). This means the form of the trial wave function in Eq. (22) is no longer applicable. In fact, consider the completeness of the coherent state, the exact wave function is a linear combination of different coherent states and should take the form

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) The ground state energy of the system with different parameters: (a) , Ω = 1, (b) , Ω = 2.

We have presented an analytical variational method for the biased quantum Rabi model. The method improves the GRWA by further considering the variational principle. The key point is that we introduce a variational parameter λ which is determined by minimizing the ground state energy. Qualitatively, the mean photon number and the correlation of the ground state are dependent on the frequency of the two-level atoms and the bias strength in our method (on the contrary, the ground state photon number and correlation in the GRWA are independent of the two-level atom and the bias, which is obviously not true). Quantitatively, comparing to the GRWA, we obtain a much more accurate ground state energy. Our approach evidently improves the GRWA in the ultra-strong coupling regime. This improvement is extremely notable if the ratio ω/Ω is small. The applicable regime of our method covers the ultra-strong coupling regime which has been achieved in the present experiment.^[9–11] Therefore, we expect that this study will be helpful for understanding the relevant experiments.