Squeezing based analytical variational method for the biased quantum Rabi model in the ultrastrong coupling regime

The quantum Rabi model (QRM) [1, 2] describes the coupling between a two-level atom and a bosonic field. Such model is of fundamental importance, since it has been widely applied in modern physics, such as in optical waveguides [3], superconducting circuits [4, 5], trapped ion [6] and so on. Given its great importance, there have been numerous theoretical studies on the Rabi model finding approximated solutions [7–9, 10–21]. Despite a great achievement of all those studies, the exact solution of QRM was only obtained recently [22–24]. This analytical solution, however, is in the form of composite transcendental function defined in power series, and it is not feasible to analyze the miscellaneous physical quantities of the QRM. Thus, theoretical research of finding the optimal analytical solution still attracts much attention.

In early days, the atom-field coupling strength of the Rabi model achieved in experiments is much smaller than the frequency of the field. Therefore, the QRM can be further simplified into the renowned Jaynes Cummings model (JCM) through rotating wave approximation (RWA), which ignores the counter-rotating-wave (CRW) terms [25]. With the tremendous development of experimental technologies, the strength of the light-matter coupling greatly increases. Recently one remarkable achievement is the realization of ultrastrong coupling [26–28]. To this end, the RWA breaks down [29, 30] since the effect of CRW terms is no longer negligible. In this respect, the generalized RWA (GRWA) [7] with a fixed displacement is proposed by considering the effect of CRW terms. However, the GRWA method does not perform well enough in the ultrastong coupling regime when the resonator frequency is small, and a more general variational method which relaxes the displacement is proposed [8, 9]. These approximation methods have been applied in the family of the QRM [8, 10–20].

The ultrastrong coupling of the QRM has been experimentally achieved in circuit quantum electrodynamics system due to operational conveniences [26–28, 30], and the so-called bias term is naturally involved. Compared to the standard form of the QRM, the additional bias term may bring in complication but more physics. For instance, the parity symmetry in the QRM is broken due to the presence of this bias. Given the interesting physics that come together with the bias term, it is therefore highly desirable to study the corresponding model analytically. Though there already exists some analytical discussions based on the GRWA and generalized variational method (GVM) for the biased QRM [9, 11], we find that there is still space to improve the accuracy in the ultrastrong coupling regime. The crucial idea to construct a close ground state in the widely used methods (e.g. GRWA, GVM) is to figure out the displacement of the oscillator. Since a quantum harmonic oscillator not only can be displaced but also can be squeezed, it is natural and necessary to consider the squeezing effect of the oscillator. The squeezing effect of the oscillator can be introduced in the same way with the displacement. The squeezing deformation of quantum state in the QRM can be naturally captured by an additional squeezing transformation applied onto the oscillator state, and thus substantial improvements may be anticipated. The squeezing plays an important role in high precision applications. For instance, it can be used to improve the sensitivity of gravity-wave detectors [31]. In quantum optics, the most common and useful way to create entangled states has been through squeezing [32–35], and it has been recognized that the intrinsic nonlinearities of atoms lead to quadrature squeezing [36–40]. With regard to this consideration, a generalized squeezing RWA (GSRWA) [15, 41] with an additional squeezing transformation was proposed to study the anisotropic QRM and shows high accuracy in large parameter regime. In this work, by taking both the displacement and squeezing transformations, we adopt the GSRWA by introducing a squeezing transformation to study the underlying physics of the biased QRM. The variational parameters in the GSRWA are determined by minimizing the ground-state energy. In order to show the improvement of our methods over the GVM, we also calculate several physical quantities by GVM and numerical method for comparison.

The paper is outlined as follows. In section 2, we give the Hamiltonian of the biased QRM and shows the details of our method. In section 3, we show the results of several physical quantities. In section 4, the applicable regime is discussed and a brief summary is finally given.

Conventionally, the Hamiltonian of the biased QRM reads [42–45]

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq_ham_0} H_B=\omega a^{\dagger}a+g\sigma_{x}(a+a^{\dagger})+\Delta\sigma_{z}+\epsilon^{\prime}\sigma_{x}, \nonumber \end{align} \tag{ 1 }$$

where $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}a^{\dagger}$ $(a)$ is the creation (annihilation) operator of the bosonic oscillator with frequency $\omega$, Pauli matrix $\sigma_{z}$ describes the two-level system with energy splitting $2\Delta$, $\newcommand{\e}{{\rm e}} \epsilon^{\prime}\sigma_{x}$ is the biased term, and g is the coupling strength between the oscillator and the qubit. Taking a rotation of the spin along y -axis and replacing the parameters with $\newcommand{\e}{{\rm e}} \epsilon^{\prime}=\epsilon/2$ and $\Delta=\Omega/2$, we can obtain

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq_ham} H&={\rm e}^{{\rm i}\frac{\pi}{4}\sigma_{y}}H_B{\rm e}^{-{\rm i}\frac{\pi}{4}\sigma_{y}} \nonumber \\ &=\omega{a}^{\dagger}a-\frac{\Omega}{2}\sigma_{x}+\frac{\epsilon}{2}\sigma_{z}+g\sigma_{z}(a^{\dagger}+a). \nonumber \end{align} \tag{ 2 }$$

For convenience, the following discussions are focusing on Hamiltonian equation (2) by default.

The model Hamiltonian in equation (2) can be divided into the main part and the perturbation as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H = H_{\rm do}+H_{\rm p}. \nonumber \end{align} \tag{ 3 }$$

Here the displaced oscillator

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{Hdo} H_{\rm do}=\omega{a}^{\dagger}a+g\sigma_{z}(a^{\dagger}+a) \nonumber \end{align} \tag{ 4 }$$

is treated as the non-perturbed Hamiltonian, and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Hpperturbation} H_{\rm p}= -\frac{\Omega}{2}\sigma_{x}+\frac{\epsilon}{2}\sigma_{z} \nonumber \end{align} \tag{ 5 }$$

can be considered as the perturbation. The ground state of the unperturbed Hamiltonian in equation (4) is doubly degenerate, which can be expressed by the direct production of a spin eigenstate and a corresponding coherent state:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle |G\rangle_{\rm do} = |\pm_z\rangle \otimes |\mp\lambda_0\rangle, \label{coh_gs0} \nonumber \end{align} \tag{ 6 }$$

where $|\pm_z\rangle$ are the eigenstates of $\sigma_z$ with eigenvalues $\pm 1$, and $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}| \pm \lambda_0 \rangle={\rm e}^{\pm \lambda_0\left(a^{\dagger}-a\right)} | 0 \rangle$ are the coherent states with $\lambda_0 = g/\omega$.

Taking the perturbation $H_{\rm p}$ into accounted, the ground state can be approximated by the linear combinations of equation (6) in the subspace of degenerate states:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{coh_gs} |G\rangle_{\rm do} = \frac{1}{\sqrt{2}}\left( |+_z\rangle \otimes |-\lambda\rangle- |-_z\rangle \otimes |\lambda\rangle\right), \nonumber \end{align} \tag{ 7 }$$

where $\lambda$ is unfixed, since $H_{\rm p}$ effectively reduces the distance of the displaced oscillator. The quantitative value of the $\lambda$ can be determined by a variational manner, namely, minimizing the ground-state energy function of the full Hamiltonian. The state in equation (7) is the starting point of the GVM, which has proven to perform better than the adiabatic approximation and the GRWA in a number of cases [8, 10, 12].

To proceed, the displacement in equation (7) can be understood through the unitary transformation on the vacuum state $|\tilde{G}\rangle_{\rm do} = U_{\rm do}$ $|G\rangle_{\rm do}$, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{udo} |\tilde{G}\rangle_{\rm do} = [\frac{1}{\sqrt{2}}\left( |+_z\rangle - |-_z\rangle \right)\otimes |0\rangle], \nonumber \end{align} \tag{ 8 }$$

where $|0\rangle$ is the vacuum Fock state, and the incorporated displacement unitary operator is

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{udo2} U_{\rm do} = {\rm e}^{-\lambda\sigma_z(a^{\dagger}-a)} . \nonumber \end{align} \tag{ 9 }$$

Correspondingly, the Hamiltonian in the transformed representation can be expressed as

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{dispH} H'=U_{\rm do}HU_{\rm do}^{\dagger}. \nonumber \end{align} \tag{ 10 }$$

Here $\lambda$ is introduced as a variational parameter rather than being fixed as $\lambda=g/\omega$ in GRWA, and this difference makes GVM obtain better accuracy. The transformed Hamiltonian reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq_ham_displace} H'=H'_0+H'_1+H'_2+H'_3, \nonumber \end{align} \tag{ 11 }$$

where

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{primeH} H'_0&= {\omega}a^{\dagger}a+\frac{\epsilon}{2}\sigma_{z}+(\omega\lambda^{2}-2\lambda{g}),\nonumber \\ H'_1&= (g-\omega\lambda)\sigma_{z}(a^{\dagger}+a),\nonumber \\ H'_2&= -\frac{\Omega}{2}\sigma_{x}{\rm cosh}[2\lambda(a^{\dagger}-a)],\nonumber \\ H'_3&= -\frac{\Omega}{2}{\rm i}\sigma_{y}{\rm sinh}[2\lambda(a^{\dagger}-a)]. \nonumber \end{align} \tag{ 12 }$$

The hyperbolic cosine and sine terms can be further expanded as [46]

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle & {\rm cosh} [2\lambda(a^{\dagger}-a)]\nonumber \\ &= F_{0}(a^{\dagger}a)+\sum_{k=0}^{\infty}[(a^{\dagger})^{2k}F_{2k+1}(a^{\dagger}a)-F_{2k}(a^{\dagger}a)a^{2k}]. \label{13} \nonumber \end{align} \tag{ 13 }$$

and

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{14} &{\rm sinh}[2\lambda(a^{\dagger}-a)]\nonumber \\ &= \sum_{k=0}^{\infty}[(a^{\dagger})^{2k+1}F_{2k+1}(a^{\dagger}a)+F_{2k+1}(a^{\dagger}a)a^{2k+1}]. \nonumber \end{align} \tag{ 14 }$$

The function F_m is defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{15} F_m(n)={\rm e}^{-2\lambda^{2}}(2\lambda)^m\frac{n!}{(n+m)!}L_{n}^{m}(4\lambda^2), \nonumber \end{align} \tag{ 15 }$$

where m and n are integers, and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{16} L_n^m(x)=\sum_{i=0}^{n}(-x)^i\frac{(n+m)!}{(m+i)!(n-i)!i!}, \nonumber \end{align} \tag{ 16 }$$

is the Laguerre polynomial. The detailed calculations see the appendix.

Similar to the displacement effect, the squeezing effect can be considered by a squeezing unitary operator

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{Us} U_{\rm s}={\rm e}^{\frac{\beta}{2}\left[a^{2}-\left(a^{\dagger}\right)^{2}\right]}. \nonumber \end{align} \tag{ 17 }$$

One can easily obtain the transformation of creation and annihilation bosonic operators as

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq_sq_op} U_{\rm s}aU_{\rm s}^{\dagger}&= a^{\dagger}{\rm sinh}\beta+a{\rm cosh}\beta,\nonumber \\ U_{\rm s}a^{\dagger}U_{\rm s}^{\dagger}&= a{\rm sinh}\beta+a^{\dagger}{\rm cosh}\beta. \nonumber \end{align} \tag{ 18 }$$

Substitute equations (18) into (11) and the transformed Hamiltonian $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\widetilde{H}=U_{\rm s}H'U_{\rm s}^{\dagger}$ is rewritten as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{dosH} \widetilde{H}=\widetilde{H}_0+\widetilde{H}_1+\widetilde{H}_2+\widetilde{H}_3, \nonumber \end{align} \tag{ 19 }$$

with

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{tildeH} \widetilde{H}_0&= \frac{\omega}{4}({\rm e}^{2\beta}+{\rm e}^{-2\beta}-2)+\frac{\omega}{2}({\rm e}^{2\beta}+{\rm e}^{-2\beta})a^{\dagger}a \nonumber \\ &\quad + \omega\lambda^2-2g\lambda+\frac{\epsilon}{2}\sigma_{z},\nonumber \\ \widetilde{H}_1&= (g-\omega\lambda)\sigma_{z}(a^\dagger+a),\nonumber \\ \widetilde{H}_2&= -\frac{\Omega}{2}\sigma_{x}{\rm cosh}[2\lambda {\rm e}^{-\beta}(a^{\dagger}-a)],\nonumber \\ \widetilde{H}_3&= -\frac{\Omega}{2}{\rm i}\sigma_{y}{\rm sinh}[2\lambda {\rm e}^{-\beta}(a^{\dagger}-a)]. \nonumber \end{align} \tag{ 20 }$$

We note that in the bosonic Fock space $\widetilde{H}_0$ only has the diagonal terms, and $\widetilde{H}_1$ shows the off-diagonal terms with the first order of a ($\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}a^{\dagger}$). By comparison, $\widetilde{H}_2$ and $\widetilde{H}_3$ are more complicated. More precisely, $\widetilde{H}_2$ can be expanded as the sum of even power of a and $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}a^{\dagger}$, while $\widetilde{H}_3$ undertakes the sum of odd power of a and $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}a^{\dagger}$. It is noteworthy that the coefficients of diagonal terms in $\widetilde{H}_0$ are evidently amplified through the squeezing transformation comparing with $H'_0$ in (12), while the off-diagonal contributions in equation (20) are either unchanged or reduced. Especially for the variational parameter $\lambda=g/\omega$ chosen in GRWA, $\widetilde{H}_1$ is vanishing. For the ultra-strong coupling regime we are interested in, we only consider diagonal terms of the Hamiltonian (19) in the following calculations, which leads to an effective adiabatic Hamiltonian:

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{eq_ham_eff} \widetilde{H}_{\rm ef\,\!f}&=\frac{\omega}{4}(\gamma+1/\gamma-2)+\frac{\omega}{2}(\gamma+1/\gamma)a^{\dagger}a\nonumber \\ &\quad +\omega\lambda^2-2g\lambda+\frac{\epsilon}{2}\sigma_{z}-\frac{\Omega}{2}\sigma_{x}F^{\prime}_0(N), \nonumber \end{align} \tag{ 21 }$$

where $F^{\prime}_0(N)={\rm e}^{-2\lambda^2{\rm e}^{-2\beta}}L^0_N(4\lambda^2{\rm e}^{-2\beta})$ and $\gamma={\rm e}^{-2\beta}$ is defined for convenience. It is obvious that the Hamiltonian in equation (21) is analytically solvable. The basis can be chosen as $|\pm_z,n\rangle$, where $|\pm_z\rangle$ are the eigenstates of $\sigma_{z}$ with eigenvalues $\pm1$, and $|n\rangle$ is the Fock state. The Hilbert space can be divided into different subspaces, in which the effective Hamiltonian equation (21) can be recast into the following form:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq_ham_matrix} \widetilde{H}_n(\lambda,\gamma)=\left( \begin{array}{@{}cc@{}} \epsilon_n^{-} & R_n \nonumber \\ R_n & \epsilon_n^{+} \end{array} \right), \nonumber \end{align} \tag{ 22 }$$

where $\newcommand{\e}{{\rm e}} \epsilon_n^{\pm}=\frac{\omega}{4}(\gamma+1/\gamma-2)+\omega(\gamma+1/\gamma)n/2+\omega\lambda^2-2g\lambda\pm \epsilon/2$, and $R_n=- \Omega F^{\prime}_0(n)/2$. From equation (22), the eigenstates and eigenvalues depending on the displacement variable $\lambda$ and the squeezing variable $\gamma$ can be easily obtained as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq_spectra} E_n^{\pm}(\lambda,\gamma)=\frac{\epsilon_n^{-}+\epsilon_n^{+}\pm\sqrt{{(\epsilon_n^{-}-\epsilon_n^{+})}^2+4R_n^2}}{2}, \nonumber \end{align} \tag{ 23 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle |\widetilde{\phi}_n^{\pm}(\lambda,\gamma)\rangle=A_n^{\pm}|-z,n\rangle+B_n^{\pm}|+z,n\rangle, \nonumber \end{align} \tag{ 24 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A_n^{\pm}&= \pm\sqrt{\frac{1}{2}\left[1\pm\frac{\epsilon_n^{-}-\epsilon_n^{+}}{\sqrt{{\left(\epsilon_n^{-}-\epsilon_n^{+}\right)}^2+4R_n^2}}\right]},\nonumber \\ B_n^{\pm}&= \sqrt{\frac{1}{2}\left[1\mp\frac{\epsilon_n^{-}-\epsilon_n^{+}}{{\left(\epsilon_n^{-}-\epsilon_n^{+}\right)}^2+4R_n^2}\right]}. \nonumber \end{align} \tag{ 25 }$$

The ground-state energy reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq_E_gs} E_{\rm G}(\lambda,\gamma)&= \frac{\omega}{4}(\gamma+1/\gamma-2)+\omega\lambda^2-2g\lambda \nonumber \\ &\quad - \frac{1}{2}\sqrt{\epsilon^2+\Omega^2{\rm e}^{-4\lambda^2\gamma}}. \nonumber \end{align} \tag{ 26 }$$

Taking $\gamma=1$, the obtained energy spectra in equation (26) recover the GVM results [9]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{\rm G}^{\rm GVM}(\lambda) = \omega \lambda^2 -2g\lambda-\sqrt{\epsilon^2+\Omega^2{\rm e}^{-4\lambda^2}} /2. \nonumber \end{align} \tag{ 27 }$$

If one further set the variational parameter $\lambda = g/\omega$, the GRWA result is recovered as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{\rm G}^{\rm GRWA} = -g^2/\omega-\frac{1}{2}\sqrt{\epsilon^2+\Omega^2{\rm e}^{-4(g/\omega)^2}} /2. \nonumber \end{align} \tag{ 28 }$$

The variational parameters $\lambda$ and $\gamma$ in the GSRWA vary with $\omega$, $\Omega$ and $\newcommand{\e}{{\rm e}} \epsilon$, and they can be determined by minimizing the $E_{\rm G}(\lambda,\gamma)$ by solving $\frac{\partial E_{\rm G}}{\partial\lambda}=0$ and $\frac{\partial E_{\rm G}}{\partial\gamma}=0$. These two nonlinear equations can not be easily solved in the analytical way. However, in the ultratrong coupling regime we interested in, the approximate solutions are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq_lambda} \lambda\approx\frac{g}{\omega+\Omega^2{\rm e}^{-4\lambda_0^2}/\sqrt{\epsilon^2+\Omega^2{\rm e}^{-4\lambda_0^2}}}, \nonumber \end{align} \tag{ 29 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq_gamma} \gamma\approx 1+\delta, \nonumber \end{align} \tag{ 30 }$$

where $\newcommand{\e}{{\rm e}} \lambda_0 = \frac{g}{\omega+\Omega^2/\sqrt{\Omega^2+\epsilon^2}}$ and $\newcommand{\e}{{\rm e}} \delta = \frac{-2\Omega^2\lambda^2}{\omega\sqrt{\epsilon^2+\Omega^2}}$. With determined values of $\lambda$ and $\gamma$, we finally obtain the analytical formula of the ground-state wavefunction,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{ground_state} {|\tilde{G}\rangle } = |\widetilde{\phi}_0^{-}(\lambda,\gamma)\rangle=-\sqrt{\frac{1}{2}[1+\frac{\epsilon}{\sqrt{\epsilon^2+ \Omega^2(F'_0(0))^2}}]}|-_z,0\rangle+ \sqrt{\frac{1}{2}[1-\frac{\epsilon}{\sqrt{\epsilon^2+\Omega^2(F'_0(0))^2}}]}|+_z,0\rangle, \nonumber \end{align} \tag{ 31 }$$

and the physical observables can thus be calculated.

Before going on, we would like to discuss the difference between our analytical method and the GVM. Equations (29) and (30) show that even if the parameter $\lambda$ takes the same value as that in the GVM, another parameter $\gamma$ may play an essential role in lowering the ground-state energy. When g and $\Omega$ are small enough, the optimized $\gamma$ in our method approaches to 1, which is the exact value presumed in the GVM. However, when g and $\Omega$ increase, the optimal value of $\gamma$ will no longer be 1. In other words, our method improves the precision of variational method when g and $\Omega$ are sizeable.

With the analytical ground-state wave function equation (31) in hand, in order to show the accuracy of the GSRWA quantitatively, we calculate several physical observables to compare with the results obtained by the GVM. The numerically exact diagonalized results are also involved as a benchmark. The quantum average of a physical observable Q is $\langle Q\rangle=\langle G|Q|G\rangle$, where $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}|G\rangle=U_{\rm do}^\dagger$ $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}U_{\rm s}^\dagger$ $|\tilde{G}\rangle$. In this work, we take the ground-state energy $E_{\rm G}=\langle H\rangle$, the mean photon number $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\langle a^{\dagger}a\rangle$, the spin orientation $\langle\sigma_{x}\rangle$, and the correlation $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\langle\sigma_{z}(a^{\dagger}+a)\rangle$ for illustration. For convenience, we set $\omega=1$ as an energy unit.

The ground-state energy obtained is given in equation (26). Figures 1(a) and (b) show the ground-state energy of the biased Rabi model as a function of g for different parameters. We see that both the results obtained by the GSRWA and the GVM agree well with the numerical ones for small coupling strength g. With the increase of g, the deviations between the GVM results and the numerical ones enlarges. In contrast, our results fit quite well with the exact ones. Figures 1(c) and (d) show the ground-state energy of the system as a function of $\Omega$ and $\newcommand{\e}{{\rm e}} \epsilon$. It is obvious that our results are better than the GVM. From these comparisons, we find that the GSRWA works much better for large g and large $\Omega$ than the GVM. The precision has been greatly improved in a wide parameter regime.

Zoom In Zoom Out Reset image size

Next we focus on the mean photon number

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \langle a^{\dagger}a\rangle=\frac{1}{4}(\gamma+1/\gamma-2)+\lambda^2, \nonumber \end{align} \tag{ 32 }$$

the spin orientation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle\sigma_{x}\rangle=\frac{\Omega {\rm e}^{-4\lambda^2\gamma}}{\sqrt{\epsilon^2+\Omega^2{\rm e}^{-4\lambda^2\gamma}}}, \nonumber \end{align} \tag{ 33 }$$

and the qubit-cavity correlation

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \langle\sigma_{z}(a^{\dagger}+a)\rangle=-2\lambda. \nonumber \end{align} \tag{ 34 }$$

As is shown in figure 2, it is obvious that our results is quite close to the numerical results for large g, where there are clear deviations for the GVM results.

Zoom In Zoom Out Reset image size

Another interesting quantity characterizing the squeezing effect is the square of the squeezing strength

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \langle (a^\dagger-a)^2 \rangle=-\gamma. \nonumber \end{align} \tag{ 35 }$$

In figure 3, we make the comparison of $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\langle (a^\dagger-a){}^2 \rangle$ calculated by different methods. One observes the result obtained by GVM is a constant  −1, but in our method, $\gamma$ changes with other parameters. The result of $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\langle (a^\dagger-a){}^2 \rangle$ through our method is in good agreement with the numerical results when $\newcommand{\e}{{\rm e}} \Omega/\epsilon$ is large, which indicates that it is necessary to take the squeezing effect into account, and our GSRWA method is a convincing way to reveal the squeezing effect of the ground-state wave function.

Zoom In Zoom Out Reset image size

As an approximation, it is necessary to discuss the applicable range of the method we proposed. It is shown that the GSRWA work quite satisfyingly for large g. According to the so-called polaron picture [47], in such ultrastrong coupling regime, the ground-state energy can be well approximated by a single oscillator. Although some numerical methods have already considered the squeezing effect [48], the analytical results are relatively rare for an intuitive picture The parameter $\lambda$ describes the displacement effect. From equation (29), we see that the optimal value of $\lambda$ is independent of squeezing rate $\gamma$, while the squeezing effect is sensitive to the rate of $\Omega/\omega$ and $\newcommand{\e}{{\rm e}} \epsilon/\omega$ shown in equation (30). When $\Omega/\omega$ is large, the $\delta$ is nonnegligible. The difference between our method and the GVM is noticeable. In such a case, our method evidently improves the GVM. When $\newcommand{\e}{{\rm e}} \epsilon/\omega$ is large, the $\delta$ is small, and it implies that a large bias reduces the squeezing.

We would like to point out that the analytical method is challenged for large $\Omega/\omega$ or large $\newcommand{\e}{{\rm e}} \epsilon/\omega$ in the sense of perturbation in equation (5). In the large $\Omega/\omega$ case, we limit our discussion to the ultrastrong coupling regime ($g \lesssim \omega$). If the coupling strength becomes much larger, one oscillator approximation (single polaron picture) is not enough. In such a case, no matter how optimal the displacement and squeezing parameters are chosen, some important physics must be missed. In such a case, a novel critical regime meets. Plenty of interesting phenomena emerge [49–52].

In conclusion, in this paper we propose an approximation method to improve the GVM to study the biased QRM by introducing the squeezing transformation. Since the QRM is an oscillator-involved model, the wave packet of its ground state is naturally expected to be compressed. In the preceding literature, various transformations to the Hamiltonian are mainly limited to the displacement effect of the wave function. The comparison of several physical quantities with different parameters show that our method works well for larger coupling strength g and atom frequency $\Omega$ than the GVM. In these parameter regimes, the effect of deformation of the oscillator state can not be neglected. Thus the additional squeezing transformation brings a prominent improvement over the GVM. A natural question arises whether such approximation method considering both the displacement and the squeezing effects works well for the low-lying excited states. We hope our study is a starting point for further investigations and has potential applications in quantum thermodynamics and quantum computing, where a very accurate approximation for obtaining energy spectra is desired [53–56].

This work is supported by Ministry of Science and Technology of China (Grant No. 2016YFA0301300), National Science Foundation of China (Grant No. 11604009, 11474211). We also acknowledge the computational support in Beijing Computational Science Research Center (CSRC).

The associated Laguerre function defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle L_{n}^{\mu}(z)=\frac{(n+\mu) !}{n ! \mu !} \sum_{l=0}^{\infty} \frac{(-n)(-n+1)(-n+2) \cdots(-n+l-1)}{l !(\mu+1)(\mu+2) \cdots(\mu+l)} z^{l}, \nonumber \end{align} \tag{ A.1 }$$

with its reduced form, namely the Laguerre function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle L_{n}(z)=L_{n}^{0}(z). \nonumber \end{align} \tag{ A.2 }$$

For the operator production, we have

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \left\{\begin{array}{@{}ll@{}}{\left(a^{\dagger}\right)^{m} a^{n}=\left(a^{\dagger}\right)^{m-n} h_{n}(\hat{N}),} & {m \geqslant n} \nonumber \\ {\left(a^{\dagger}\right)^{m} a^{n}=h_{m}(\hat{N}) a^{n-m},} & {m<n} \end{array}\right. \nonumber \end{align} \tag{ A.3 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_{n}(\hat{N})=\hat{N}(\hat{N}-1)(\hat{N}-2) \cdots(\hat{N}-n+1). \nonumber \end{align} \tag{ A.4 }$$

With the above equations in hand, we can obtain

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \cosh \nu\left(a^{\dagger}-a\right) &=\frac{1}{2}\left[{\rm e}^{\nu}\left(a^{\dagger}-a\right)+{\rm e}^{-\nu}\left(a^{\dagger}-a\right)\right] \nonumber \\ &=\frac{1}{2} {\rm e}^{-\nu^{2} / 2}\left[{\rm e}^{\nu a^{\dagger}} {\rm e}^{-\nu a}+{\rm e}^{-\nu a^{\dagger}} {\rm e}^{\nu a}\right] \nonumber \\ &=\frac{1}{2} {\rm e}^{-\nu^{2} / 2} \sum_{m, n}^{\infty} \frac{1}{m ! n !}\left[\nu^{m}(-\nu)^{n}+(-\nu)^{m} \nu^{n}\right]\left(a^{\dagger}\right)^{m} a^{n}. \nonumber \end{align} \tag{ A.5 }$$

Here, we use the formula

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm e}^{(A+B)}={\rm e}^{A} {\rm e}^{B} {\rm e}^{-\frac{1}{2}[A, B]}. \nonumber \end{align} \tag{ A.6 }$$

For $m-n=2 k \geqslant 0$, there is

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle I_{x}^{+} & =\frac{1}{2} {\rm e}^{-\nu^{2} / 2} \sum_{m, n}^{\infty} \frac{1}{m ! n !}\left[\nu^{m}(-\nu)^{n}+(-\nu)^{m} \nu^{n}\right]\left(a^{\dagger}\right)^{m} a^{n} \nonumber \\ &=\frac{1}{2} {\rm e}^{-\nu^{2} / 2} \sum_{k}^{\infty} \sum_{n}^{\infty} \frac{1}{(n+2 k) ! n !}\left[\nu^{n+2 k}(-\nu)^{n}+(-\nu)^{(n+2 k)} \nu^{n}\right]\left(a^{\dagger}\right)^{(n+2 k)} a^{n} \nonumber \\ &=\frac{1}{2} {\rm e}^{-\nu^{2} / 2} \sum_{k}^{\infty} \nu^{2 k}\left(a^{\dagger}\right)^{2 k} \sum_{n}^{m} \frac{(-)^{n} h_{n}(\hat{N})}{(n+2 k) ! n !}\left(2 \nu^{2 n}\right) \nonumber \\ &={\rm e}^{-\nu^{2} / 2} \sum_{k}^{\infty} \nu^{2 k}\left(a^{\dagger}\right)^{2 k} \frac{(\hat{N}+2 k) !}{N !} \frac{\hat{N} !}{(\hat{N}+2 k) !} \sum_{n}^{\infty} \frac{(-)^{n} h_{n}(\hat{N})}{(n+2 k) ! n !} \nu^{2 n} \nonumber \\ &={\rm e}^{-\nu^{2} / 2} \sum_{k}^{\infty} \nu^{2 k}\left(a^{\dagger}\right) \frac{(\hat{N}+2 k) !}{\hat{N} !} L_{\hat{N}}^{2 k}\left(\nu^{2}\right). \nonumber \end{align} \tag{ A.7 }$$

For $m-n= -2 k \leqslant 0$, there is

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle I_{x}^{-} &=\frac{1}{2} {\rm e}^{-\nu^{2} / 2} \sum_{m, n}^{\infty} \frac{1}{m ! n !}\left[\nu^{m}(-\nu)^{n}+(-\nu)^{m} \nu^{n}\right]\left(a^{\dagger}\right)^{m} a^{n} \nonumber \\ &={\rm e}^{-\nu^{2} / 2} \sum_{k}^{\infty} \nu^{2 k} \frac{(\hat{N}+2 k) !}{\hat{N !}} L_{N}^{2 k}\left(\nu^{2}\right) a^{2 k}. \nonumber \end{align} \tag{ A.8 }$$

With the above all and the definition of F_m(n) in equation (15), the equation (13) can be obtained after some algebra. In the similar way, equation (14) can be obtained as well.