Abstract An analytical method is developed to study the two-mode quantum Rabi model. For certain specific parameter conditions, especially for the resonant conditions, we obtain an infinite number of the exact solutions of the eigenfunctions and associated energies. It is shown that there exist new types of the exact energies which do not correspond to the level-crossings. Our analytical method may find applications in some related models. Keywords:exact solutions;two-mode quantum Rabi model;level-crossings
Since the exact solutions to the quantum Rabi model (QRM) have been presented [1–4], the search of the analytical solutions has attracted renewed interest of research in some related physical models, such as the anisotropic QRM [5, 6], the two-photon QRM [7–11], the two-mode QRM [12–16], and several related models [17–21]. These models can provide a good theoretical description of the interaction between matter and light in different situations. Parts of these models have been realized in a variety of physical systems [22, 23].
Recent studies have shown that these models have two interesting features. Firstly, for some specially chosen parameter conditions, there are an infinite number of the exact energies and associated eigenfunctions. The exact eigenfunctions are usually given in terms of elementary functions [1–3, 5]. In the two-photon QRM, certain exact eigenfunctions are even expressed by some special functions [11]. For general parameter conditions, the energies and associated eigenfunctions cannot be obtained algebraically and explicitly. Secondly, it was shown that if the underlying parity symmetries in these models are not broken by some additional terms, all the known exact energies are related to the level-crossings.
The two-mode QRM is an important extension of the QRM. In Bargmann Hilbert spaces of entire functions, the Bethe ansatz method and the algebraic method have been applied to obtain the exact parts of the energy spectrum associated with polynomial eigenfunctions [12, 15]. Later, these exact energies have been obtained by using extended squeezed states. The resulting eigenfunctions are given by a finite expansion of the extended squeezed states [13, 16].
In the present work, we develop a different analytical method to study the solutions of the two-mode QRM. The eigenfunctions are expanded as a series of the confluent hypergeometric functions (CHFs). It is shown that when the parameters satisfy certain specific conditions, the series expansion become finite, thereby resulting in an infinite number of the exact solutions of the energies. In particular, it is found that at the resonant conditions, there are also an infinite number of the exact energies. In a comparison with the previously existing exact energies, it is further found that the obtained exact results contain the previous results as well as the new results. These new exact energies are associated the non-polynomial eigenfunctions, and do not correspond to the level-crossings. Therefore, our analytical results show that there are two types of exact eigenfunctions, polynomial and non-polynomial eigenfunctions. The polynomial eigenfunctions are associated with the level-crossings, and the non-polynomial eigenfunctions are not. In addition, it is shown that the analytical method may be applied to present the exact energies and eigenfunctions of some related physical models.
2. Exact solutions of the two-mode QRM
We consider the two-mode QRM described by the following Hamiltonian (ℏ=1)$\begin{eqnarray}H=\omega ({a}_{1}^{\dagger }{a}_{1}+{a}_{2}^{\dagger }{a}_{2})+{\rm{\Delta }}{\sigma }_{z}+g{\sigma }_{x}({a}_{1}^{\dagger }{a}_{2}^{\dagger }+{a}_{1}{a}_{2}),\end{eqnarray}$where ${a}_{1,2}^{\dagger }$ and a1,2 are the creation and annihilation operators for two bosonic modes with the same frequency ω, σx,z are the usual Pauli matrices, 2Δ is the energy separation between the two levels, and g denotes the interaction strength. In the two-mode QRM, the eigenstate $| \psi \rangle $ is expanded as$\begin{eqnarray}| \psi \rangle ={\psi }_{1}({a}_{1}^{\dagger },{a}_{2}^{\dagger })| 0\rangle | \uparrow \rangle +{\psi }_{2}({a}_{1}^{\dagger },{a}_{2}^{\dagger })| 0\rangle | \downarrow \rangle .\end{eqnarray}$Here as in the QRM, ψ1,2 must be analytical entire functions of the creation operators ${a}_{1,2}^{\dagger }$. $| 0\rangle \equiv | 0{\rangle }_{1}| 0{\rangle }_{2}$ is the vacuum state for the two bosonic modes. $| \uparrow \rangle $ and $| \downarrow \rangle $ are the eigenstates of σz with eigenvalues 1 and −1. Substituting this expansion into the Schrödinger equation $H| \psi \rangle =E| \psi \rangle $, we get the eigenvalue equation$\begin{eqnarray}\begin{array}{l}([H,{\psi }_{1}]+{\psi }_{1}H-E{\psi }_{1})| 0\rangle | \uparrow \rangle \\ \quad +\,([H,{\psi }_{2}]+{\psi }_{2}H-E{\psi }_{2})| 0\rangle | \downarrow \rangle =0.\end{array}\end{eqnarray}$By applying the well-known relations, ${a}_{\mathrm{1,2}}| 0\rangle =0$, $[{a}_{1}^{\dagger },{\psi }_{\mathrm{1,2}}]=0$, $[{a}_{2}^{\dagger },{\psi }_{\mathrm{1,2}}]=0$, $[{a}_{1},{\psi }_{\mathrm{1,2}}]=\partial {\psi }_{\mathrm{1,2}}/\partial {a}_{1}^{\dagger }$, and $[{a}_{2},{\psi }_{\mathrm{1,2}}]=\partial {\psi }_{\mathrm{1,2}}/\partial {a}_{2}^{\dagger }$, together with the linear combinations, ${\varphi }_{1}=({\psi }_{1}+{\psi }_{2})/2$ and ${\varphi }_{2}=({\psi }_{1}-{\psi }_{2})/2$, we have two coupled operator-type differential equations$\begin{eqnarray}\begin{array}{l}{z}_{1}\displaystyle \frac{\partial {\varphi }_{1}}{\partial {z}_{1}}+{z}_{2}\displaystyle \frac{\partial {\varphi }_{1}}{\partial {z}_{2}}+g\displaystyle \frac{{\partial }^{2}{\varphi }_{1}}{\partial {z}_{1}\partial {z}_{2}}\\ \quad +\,({{gz}}_{1}{z}_{2}-E){\varphi }_{1}+{\rm{\Delta }}{\varphi }_{2}=0,\end{array}\end{eqnarray}$$\begin{eqnarray}\begin{array}{l}{z}_{1}\displaystyle \frac{\partial {\varphi }_{2}}{\partial {z}_{1}}+{z}_{2}\displaystyle \frac{\partial {\varphi }_{2}}{\partial {z}_{2}}-g\displaystyle \frac{{\partial }^{2}{\varphi }_{2}}{\partial {z}_{1}\partial {z}_{2}}\\ \quad -({{gz}}_{1}{z}_{2}+E){\varphi }_{2}+{\rm{\Delta }}{\varphi }_{1}=0,\end{array}\end{eqnarray}$where ${z}_{1}={a}_{1}^{\dagger }$, ${z}_{2}={a}_{2}^{\dagger }$, and we have set ω=1 for brevity. It is evident that all the terms in this equation are mutually commutable. Aa a result, one can formally regard them as c-number differential equations [24, 25]. In addition, the two-mode QRM has a conserved operator$\begin{eqnarray}C=-\displaystyle \frac{1}{4}({n}_{1}-{n}_{2}+1)({n}_{1}-{n}_{2}-1),\end{eqnarray}$where ${n}_{1}={a}_{1}^{\dagger }{a}_{1}$ and ${a}_{2}^{\dagger }{a}_{2}$ are the photon-number operators of mode 1 and mode 2, respectively. This conserved operator allows us to decouple the Fock–Hilbert space into the direct sum of infinite number of subspaces labeled by $\kappa =({n}_{1}-{n}_{2}+1)/2$ [12, 15]. Here κ is usually called as the Bargmann index, and can take any positive integers or half-integers, $\kappa =1/2,1,3/2,\ldots $. Due to the presence of the conserved operator C, φ1,2 can take the form$\begin{eqnarray}{\varphi }_{\mathrm{1,2}}={z}_{1}^{2\kappa -1}{f}_{\mathrm{1,2}}(z),\end{eqnarray}$with $z={z}_{1}{z}_{2}$. After a simple calculation, we have the coupled equations for ${f}_{\mathrm{1,2}}(z)$$\begin{eqnarray}\begin{array}{l}{gz}\displaystyle \frac{{{\rm{d}}}^{2}{f}_{1}}{{\rm{d}}{z}^{2}}+(2\kappa g+2z)\displaystyle \frac{{\rm{d}}{f}_{1}}{{\rm{d}}z}\\ \quad +\,(2\kappa -1-E+{gz}){f}_{1}+{\rm{\Delta }}{f}_{2}=0,\end{array}\end{eqnarray}$$\begin{eqnarray}\begin{array}{l}{gz}\displaystyle \frac{{{\rm{d}}}^{2}{f}_{2}}{{\rm{d}}{z}^{2}}+(2\kappa g-2z)\displaystyle \frac{{\rm{d}}{f}_{2}}{{\rm{d}}z}\\ \quad +\,(1-2\kappa +E+{gz}){f}_{2}-{\rm{\Delta }}{f}_{1}=0.\end{array}\end{eqnarray}$It is easily found that if (f1(z), f2(z)) is a solution, $({f}_{2}(-z),{f}_{1}(-z))$ is also a solution. This property can simplify our discussions. If we further take the following transformations$\begin{eqnarray}{f}_{\mathrm{1,2}}(z)={{\rm{e}}}^{-\eta z}{\phi }_{\mathrm{1,2}}(y),y=\xi z,\end{eqnarray}$we get$\begin{eqnarray}y\displaystyle \frac{{{\rm{d}}}^{2}{\phi }_{1}}{{\rm{d}}{y}^{2}}+(2\kappa +{\lambda }_{1}^{+}y)\displaystyle \frac{{\rm{d}}{\phi }_{1}}{{\rm{d}}y}+({\lambda }_{2}^{+}+{\lambda }_{3}^{+}y){\phi }_{1}+\displaystyle \frac{{\rm{\Delta }}}{g\xi }{\phi }_{2}=0,\end{eqnarray}$$\begin{eqnarray}y\displaystyle \frac{{{\rm{d}}}^{2}{\phi }_{2}}{{\rm{d}}{y}^{2}}+(2\kappa +{\lambda }_{1}^{-}y)\displaystyle \frac{{\rm{d}}{\phi }_{2}}{{\rm{d}}y}+({\lambda }_{2}^{-}+{\lambda }_{3}^{-}y){\phi }_{2}-\displaystyle \frac{{\rm{\Delta }}}{g\xi }{\phi }_{1}=0,\end{eqnarray}$where$\begin{eqnarray}{\lambda }_{1}^{\pm }=\displaystyle \frac{-2\eta g\pm 2}{g\xi },\end{eqnarray}$$\begin{eqnarray}{\lambda }_{2}^{\pm }=\displaystyle \frac{-2\eta g\kappa \pm (2\kappa -1-E)}{g\xi },\end{eqnarray}$$\begin{eqnarray}{\lambda }_{3}^{\pm }=\displaystyle \frac{g{\eta }^{2}\mp 2\eta +g}{g{\xi }^{2}}.\end{eqnarray}$Here the two parameters η and ξ are constants to be determined.
The power series expansion is a conventional method to solve equations (11) and (12) by expanding φ1,2 by power series of y with unknown coefficients. This method has been used for the QRM and its related models. For a general parameter condition, one has a higher-order recurrence relation for the coefficients. If the model parameters satisfy certain special condition, the recurrence relation can be terminated. Such termination leads to the exact energies in many Rabi-related models. In the two-mode QRM, by using the power series expansion, we get the exact energies [12, 15]$\begin{eqnarray}\displaystyle \frac{{E}_{N,\kappa }^{{\rm{I}}}}{\omega }=-1+2(N+\kappa )\sqrt{1-\displaystyle \frac{{g}^{2}}{{\omega }^{2}}}.\end{eqnarray}$Here N≥1 , ${\eta }_{\pm }=\pm (\omega -\sqrt{{\omega }^{2}-{g}^{2}})/g$, and the parameter ξ can be chosen as an arbitrary non-zero constant. This exact energies were first found in [12, 15] via different methods, and correspond to the level-crossings in the energy spectrum [13]. The exact energies exist for certain specific parameter relations. For example, in the cases of N=1, 2, the parameter relations are given as$\begin{eqnarray}\displaystyle \frac{{{\rm{\Delta }}}^{2}}{{\omega }^{2}}+4(2\kappa +1)\displaystyle \frac{{g}^{2}}{{\omega }^{2}}-4=0,\end{eqnarray}$$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{{{\rm{\Delta }}}^{4}}{{\omega }^{4}}+(4(7+6\kappa )\displaystyle \frac{{g}^{2}}{{\omega }^{2}}-20)\displaystyle \frac{{{\rm{\Delta }}}^{2}}{{\omega }^{2}}+128{\left(1+\kappa )\displaystyle \frac{{g}^{2}}{{\omega }^{2}}-1\right)}^{2}\\ \quad +\,64((1+\kappa )\displaystyle \frac{{g}^{4}}{{\omega }^{4}}-1)=0.\end{array}\end{eqnarray}$It is evident that if we plot the energies as a function of g/ω with the fixed values of Δ/ω and κ, the exact energies (16) appear in some isolated values of g/ω [13].
In this work, our aim is to find new expansion functions which allow us to analytically obtain more exact energies. It is observed that under the condition of ${\lambda }_{3}^{\pm }=0$ and Δ=0, equations (11) and (12) decouple, and have a similar form with the confluent hypergeometric equation [26]. This motivates us to search for solutions with a series expansion in terms of the CHFs. The series expansion in terms of the CHFs has been used in different systems [27–30]. We shall show that under certain special parameter conditions, the series expansion can be terminated. This gives the exact analytical results of the energies. In the following, we discuss the series expansion of the solutions of equations (11) and (12) in terms of the CHFs.
2.1. Series expansion in terms of CHFs
The series expansion of the solutions to equations (11) and (12) in terms of the CHFs is given as$\begin{eqnarray}{\phi }_{1}(y)=\sum _{n=0}^{\infty }{a}_{n}\ F({\alpha }_{n},\gamma ,y),\ \ {\phi }_{2}(y)=\sum _{n=0}^{\infty }{b}_{n}F({\alpha }_{n},\gamma ,y),\end{eqnarray}$where F(αn, γ, y) is the CHF with ${\alpha }_{n}={\alpha }_{0}+n$. Here α0 and γ are constants to be determined. After substituting these expansions into equations (11) and (12), and collecting the like terms, we obtain the recurrence relation for the coefficients an and bn (see details in appendix A)$\begin{eqnarray}{A}_{n+1}{W}_{n+1}+{B}_{n}{W}_{n}+{C}_{n-1}{W}_{n-1}=0,\end{eqnarray}$with$\begin{eqnarray}{W}_{n}=\left(\begin{array}{c}{a}_{n}\\ {b}_{n}\end{array}\right).\end{eqnarray}$It is a hard task to obtain an analytical solution for this higher-order recurrence relation. However under the condition$\begin{eqnarray}{W}_{n}=0\ \ \mathrm{if}\ \ n\gt N\end{eqnarray}$one may cut off the recurrence expressions, and solve them step by step. Such a termination condition may allow us to find the exact analytical results of the energies.
2.2. Exact energies
From the terminated condition (22), we have the explicit expressions of the exact energies (see details in appendix B)$\begin{eqnarray}\displaystyle \frac{{E}_{N,\kappa }^{\mathrm{II}}}{\omega }=-1+N\sqrt{1-\displaystyle \frac{{g}^{2}}{{\omega }^{2}}}.\end{eqnarray}$Here N≥1, ${\eta }_{\pm }={\xi }_{\pm }/2=\pm (\omega -\sqrt{{\omega }^{2}-{g}^{2}})/g$, ${\alpha }_{0}\,=\kappa -N/2$, and γ=2κ. The exact energies exist for certain specific parameter relations. For example, in the cases of N=1, 2, 3, we have$\begin{eqnarray}\left(\displaystyle \frac{{\rm{\Delta }}}{\omega }-2\kappa +1\right)\left(\displaystyle \frac{{\rm{\Delta }}}{\omega }+2\kappa -1\right)=0,\end{eqnarray}$$\begin{eqnarray}\displaystyle \frac{{{\rm{\Delta }}}^{4}}{{\omega }^{4}}+4\left(\displaystyle \frac{{g}^{2}}{{\omega }^{2}}-1+2(1-\kappa )\kappa \right)\displaystyle \frac{{{\rm{\Delta }}}^{2}}{{\omega }^{2}}+16{\left(1-\kappa \right)}^{2}{\kappa }^{2}=0,\end{eqnarray}$$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{{{\rm{\Delta }}}^{3}}{{\omega }^{3}}\pm \left(1-2\kappa \right)\displaystyle \frac{{{\rm{\Delta }}}^{2}}{{\omega }^{2}}+(8\displaystyle \frac{{g}^{2}}{{\omega }^{2}}-4-{\left(1-2\kappa \right)}^{2})\displaystyle \frac{{\rm{\Delta }}}{\omega }\\ \quad \pm \,(4{\kappa }^{2}-1)(2\kappa -3)=0.\end{array}\end{eqnarray}$
Now we discuss the relation between the two exact energies ${E}_{N,\kappa }^{{\rm{I}}}$ and ${E}_{N,\kappa }^{\mathrm{II}}$. Due to $2(N+\kappa )\geqslant 3$ in ${E}_{N,\kappa }^{{\rm{I}}}$, it follows that the exact energies ${E}_{N,\kappa }^{\mathrm{II}}$ with N=1,2 are beyond the exact energies ${E}_{N,\kappa }^{{\rm{I}}}$. In addition, for certain given values of N and κ, the resulting CHFs become polynomial. In such situations, we get$\begin{eqnarray}{E}_{N,\kappa }^{{\rm{I}}}={E}_{2(N+\kappa ),\kappa }^{\mathrm{II}}.\end{eqnarray}$This result indicates that all the exact energies ${E}_{N,\kappa }^{{\rm{I}}}$ are contained in the exact energies ${E}_{N,\kappa }^{\mathrm{II}}$. For example, from the two parameter relations (17) and (26), we have ${E}_{1,1/2}^{{\rm{I}}}={E}_{3,1/2}^{\mathrm{II}}$. But the exact energies ${E}_{3,\kappa }^{\mathrm{II}}$ with $\kappa \gt 1/2$ are not included in the exact energies ${E}_{N,\kappa }^{{\rm{I}}}$. These obtained new exact energies are not found previously. In the following, we shall show that they do not correspond to the level-crossings.
We first study the interesting exact energies ${E}_{N,\kappa }^{\mathrm{II}}$ with N=1, since the parameter g/ω does not appear in the parameter relation (24). Without loss of generality, we take Δ>0. From the parameter relation (24), we have the parameter relation Δ/ω=(2κ−1) corresponding to the resonant conditions. This means that at the resonant conditions, there exist the exact energies ${E}_{1,\kappa }^{\mathrm{II}}/\omega =-1+\sqrt{1-{g}^{2}/{\omega }^{2}}$. The corresponding solutions are not expressed in terms of elementary functions (see details in appendix C). We note that in the two-photon QRM, there also exists the exact energy at the resonant condition (see details in appendix D). The reason behind this deserves to be studied in future. In figure 1, we show the numerical and exact analytical results of the energy spectrum of the lowest energy levels as a function of g/ω with ω=1 and Δ/ω=1. This choice of Δ/ω=1 corresponds to the resonant condition of Δ/ω=(2κ−1) with κ=1. Here we discuss four sets of the exact energies, marked by square, triangles, circles, and dots. The square is for the exact energy ${E}_{2,1}^{\mathrm{II}}$. The triangles denote the exact energies ${E}_{3,\kappa }^{\mathrm{II}}$ with $\kappa =1/2,1,3/2$. The circles denote the exact energies ${E}_{1,1}^{\mathrm{II}}$. The dots denote the exact energies ${E}_{1,\kappa }^{{\rm{I}}}$ with κ=1/2, 1. It is seen from the analytical and numerical results that the exact energies ${E}_{2,1}^{\mathrm{II}}$ and ${E}_{3,3/2}^{\mathrm{II}}$ do not correspond to the level-crossings. Furthermore, ${E}_{1,1}^{\mathrm{II}}$ shows a continuous dependence on the parameter g/ω up to g/ω=1 where it meets ${E}_{3,1}^{\mathrm{II}}$. In addition, as shown figure 2, we observe that depending on the chosen values of Δ/ω and κ, the position of the exact energies ${E}_{1,\kappa }^{\mathrm{II}}$ changes. Therefore, these numerical and analytical results demonstrate that the two-mode QRM supports new exact energies which are not associated with the level-crossings.
Figure 1.
New window|Download| PPT slide Figure 1.Energy spectrum of the two-mode QRM as a function of g/ω with ω=1 and Δ/ω=1. The square denotes the exact energy ${E}_{2,1}^{\mathrm{II}}$. The triangles denote the exact energies ${E}_{3,\kappa }^{\mathrm{II}}$ with $\kappa =1/2,1,3/2$. The circles denote the exact energies ${E}_{1,1}^{\mathrm{II}}$. The dots denote the exact energies ${E}_{1,\kappa }^{{\rm{I}}}$ with $\kappa =1/2,1$.
Figure 2.
New window|Download| PPT slide Figure 2.Energy spectrum of the two-mode QRM as a function of g/ω for two different values of (a) Δ/ω=2 and (b) Δ/ω=3 with ω=1. In (a), the circles denote the exact energies ${E}_{1,3/2}^{\mathrm{II}}$ at the resonant condition of ${\rm{\Delta }}/\omega =2\kappa -1=2$, and in (b), the circles denote the exact energies ${E}_{1,2}^{\mathrm{II}}$ at the resonant condition of ${\rm{\Delta }}/\omega =2\kappa -1=3$.
3. Conclusion
In conclusion, the series expansion of the eigenfunctions in terms of the CHFs has been introduced to study the two-mode QRM. The conditions for the finite series expansion have been applied to find an infinite number of the exact energies. Depending on whether they are associated with the level-crossings, two different types of the exact energies are identified. They correspond to the polynomial and non-polynomial eigenfunctions, respectively. In particular, it has been found that there exist the exact energies at the resonant conditions.
This analytical method has been applied to present the exact energies for the two-photon QRM. It is believed that our analytical method could be used to find new exact energies associated with the non-polynomial eigenfunctions in some related physical models, such as the anisotropic two-photon QRM, and the anisotropic two-mode QRM [16]. In addition, one could use the new exact energies to discuss the spectral collapse phenomenon [13]. It is an interesting problem to explore whether there is a hidden Lie algebraic structure which explains the existence of the non-polynomial eigenfunctions in the two-mode QRM [15].
Acknowledgments
Supported by the National Natural Science Foundation of China under Grant No. 11965011.
Appendix A. Derivation of the recurrence relation (20)
In this section, we present a detailed derivation of the recurrence relation (20). Substituting the series expansion ${\phi }_{1}(y)={\sum }_{n=0}^{\infty }{a}_{n}{u}_{n}(y)$ and ${\phi }_{2}(y)={\sum }_{n=0}^{\infty }{b}_{n}{u}_{n}(y)$ with ${u}_{n}(y)=F({\alpha }_{n},\gamma ,y)$ into equations (11) and (12) gives$\begin{eqnarray}\displaystyle \begin{array}{l}\sum _{n=0}^{\infty }{a}_{n}\left(y\displaystyle \frac{{{\rm{d}}}^{2}{u}_{n}}{{\rm{d}}{y}^{2}}+(2\kappa +{\lambda }_{1}^{+}y)\displaystyle \frac{{\rm{d}}{u}_{n}}{{\rm{d}}y}\right.\\ \quad \left.+\,({\lambda }_{2}^{+}+{\lambda }_{3}^{+}y){u}_{n}\right)+\displaystyle \frac{{\rm{\Delta }}}{g\xi }\sum _{n=0}^{\infty }{b}_{n}{u}_{n}=0,\end{array}\end{eqnarray}$$\begin{eqnarray}\displaystyle \begin{array}{l}\sum _{n=0}^{\infty }{b}_{n}\left(y\displaystyle \frac{{{\rm{d}}}^{2}{u}_{n}}{{\rm{d}}{y}^{2}}+(2\kappa +{\lambda }_{1}^{-}y)\displaystyle \frac{{\rm{d}}{u}_{n}}{{\rm{d}}y}\right.\\ \quad +\,\left.({\lambda }_{2}^{-}+{\lambda }_{3}^{-}y){u}_{n}\right)-\displaystyle \frac{{\rm{\Delta }}}{g\xi }\sum _{n=0}^{\infty }{a}_{n}{u}_{n}=0.\end{array}\end{eqnarray}$We use the following relations [26]$\begin{eqnarray}y\displaystyle \frac{{{\rm{d}}}^{2}{u}_{n}}{{\rm{d}}{y}^{2}}=(y-\gamma )\displaystyle \frac{{\rm{d}}{u}_{n}}{{\rm{d}}y}+{\alpha }_{n}{u}_{n},\end{eqnarray}$$\begin{eqnarray}y\displaystyle \frac{{\rm{d}}{u}_{n}}{{\rm{d}}y}={\alpha }_{n}({u}_{n+1}-{u}_{n}),\end{eqnarray}$$\begin{eqnarray}{{yu}}_{n}=({\alpha }_{n}-\gamma ){u}_{n-1}+(\gamma -2{\alpha }_{n}){u}_{n}+{\alpha }_{n}{u}_{n+1}.\end{eqnarray}$To further simplify the above equations$\begin{eqnarray}\displaystyle \begin{array}{l}\sum _{n=0}^{\infty }{a}_{n}\left(2\kappa -\gamma )\displaystyle \frac{{\rm{d}}{u}_{n}}{{\rm{d}}y}+{Q}_{n}^{+}{u}_{n+1}+{P}_{n}^{+}{u}_{n}+{R}_{n}^{+}{u}_{n-1}\right)\\ \quad +\,\displaystyle \frac{{\rm{\Delta }}}{g\xi }\sum _{n=0}^{\infty }{b}_{n}{u}_{n}=0,\end{array}\end{eqnarray}$$\begin{eqnarray}\displaystyle \begin{array}{l}\sum _{n=0}^{\infty }{b}_{n}\left(2\kappa -\gamma )\displaystyle \frac{{\rm{d}}{u}_{n}}{{\rm{d}}y}+{Q}_{n}^{-}{u}_{n+1}+{P}_{n}^{-}{u}_{n}+{R}_{n}^{-}{u}_{n-1}\right)\\ \quad -\,\displaystyle \frac{{\rm{\Delta }}}{g\xi }\sum _{n=0}^{\infty }{a}_{n}{u}_{n}=0,\end{array}\end{eqnarray}$with$\begin{eqnarray*}\begin{array}{rcl}{R}_{n}^{\pm } & = & {\lambda }_{3}^{\pm }({\alpha }_{n}-\gamma ),\\ {P}_{n}^{\pm } & = & {\lambda }_{3}^{\pm }(\gamma -2{\alpha }_{n})+{\lambda }_{2}^{\pm }-{\lambda }_{1}^{\pm }{\alpha }_{n},\\ {Q}_{n}^{\pm } & = & (1+{\lambda }_{1}^{\pm }+{\lambda }_{3}^{\pm }){\alpha }_{n}.\end{array}\end{eqnarray*}$Since ${\rm{d}}{u}_{n}/{\rm{d}}y$ cannot be expressed as a linear combination of functions un, we demand $\gamma =2\kappa $. After equating coefficient of un(y)=F(αn, γ, y), we have$\begin{eqnarray}{R}_{n+1}^{+}{a}_{n+1}+{P}_{n}^{+}{a}_{n}+\displaystyle \frac{{\rm{\Delta }}}{g\xi }{b}_{n}+{Q}_{n-1}^{+}{a}_{n-1}=0,\end{eqnarray}$$\begin{eqnarray}{R}_{n+1}^{-}{b}_{n+1}+{P}_{n}^{-}{b}_{n}-\displaystyle \frac{{\rm{\Delta }}}{g\xi }{a}_{n}+{Q}_{n-1}^{-}{b}_{n-1}=0.\end{eqnarray}$If we let$\begin{eqnarray}{W}_{n}=\left(\begin{array}{c}{a}_{n}\\ {b}_{n}\end{array}\right),\end{eqnarray}$we have the recurrence relation (20)$\begin{eqnarray}{A}_{n+1}{W}_{n+1}+{B}_{n}{W}_{n}+{C}_{n-1}{W}_{n-1}=0,\end{eqnarray}$where An, Bn and Cn are given as$\begin{eqnarray*}{A}_{n}=\left(\begin{array}{cc}{R}_{n}^{+} & 0\\ 0 & {R}_{n}^{-}\end{array}\right),{B}_{n}=\left(\begin{array}{cc}{P}_{n}^{+} & \tfrac{{\rm{\Delta }}}{g\xi }\\ -\tfrac{{\rm{\Delta }}}{g\xi } & {P}_{n}^{-}\end{array}\right),{C}_{n}=\left(\begin{array}{cc}{Q}_{n}^{+} & 0\\ 0 & {Q}_{n}^{-}\end{array}\right).\end{eqnarray*}$
Appendix B. Derivation of the exact energies in equation (23)
In this section, a detailed derivation of the exact energies (23) is given. By setting n=N+1 in the recurrence relation (20), we have ${C}_{N}{W}_{N}=0$. For non-zero solution of aN and bN, the determinant of the matrix CN must vanish, thereby resulting in the condition, $(1+{\lambda }_{1}^{+}+{\lambda }_{3}^{+})(1+{\lambda }_{1}^{-}+{\lambda }_{3}^{-}){\alpha }_{N}^{2}=0$. If we choose ${\alpha }_{N}=0$. we get ${\alpha }_{0}=-N$ corresponding to the power series expansion. We do not consider this situation. By setting n=−1 in the recurrence relation (20), we have ${A}_{0}{W}_{0}=0$, and thus obtain the condition, ${\lambda }_{3}^{+}{\lambda }_{3}^{-}{({\alpha }_{0}-\gamma )}^{2}=0$. After a detailed calculation, it is found that for a physically acceptable solutions, we take two sets of conditions, (1) ${\lambda }_{3}^{-}=0$, $1+{\lambda }_{1}^{+}+{\lambda }_{3}^{+}=0$, ${a}_{0}=0$, and bN=0, and (2) ${\lambda }_{3}^{+}=0$, $1+{\lambda }_{1}^{-}+{\lambda }_{3}^{-}=0$, aN=0, and ${b}_{0}=0$. For the first set of conditions, by setting n=0 and n=N in the recurrence relation (20), we further get ${\lambda }_{2}^{-}-{\lambda }_{1}^{-}{\alpha }_{0}=0$ and ${\lambda }_{3}^{+}(\gamma -2{\alpha }_{N})+{\alpha }_{N}+{\lambda }_{2}^{+}-(1+{\lambda }_{1}^{+}){\alpha }_{N}\,=0$. For the second set of conditions, by setting n=0 and n=N in the recurrence relation (20), we have ${\lambda }_{2}^{+}-{\lambda }_{1}^{+}{\alpha }_{0}=0$ and ${\lambda }_{3}^{+}(\gamma -2{\alpha }_{N})+{\lambda }_{2}^{+}-{\lambda }_{1}^{+}{\alpha }_{N}=0$. The parameter η is determined by the conditions, ${\lambda }_{3}^{\pm }=0$. After substituting η into the conditions, $1+{\lambda }_{1}^{\mp }+{\lambda }_{3}^{\mp }=0$, the parameters ξ are given. Once the parameters η and ξ are determined, the parameters ${\alpha }_{0}$ and the energies E are given by the conditions, ${\lambda }_{2}^{\mp }-{\lambda }_{1}^{\mp }{\alpha }_{0}=0$ and ${\lambda }_{3}^{\pm }(\gamma -2{\alpha }_{N})+{\lambda }_{2}^{\pm }-{\lambda }_{1}^{\pm }{\alpha }_{N}=0$.
In the first situation with ${\lambda }_{3}^{-}=0$, we have the parameters and the energies$\begin{eqnarray}\begin{array}{rcl}{\eta }_{-} & = & \displaystyle \frac{{\xi }_{-}}{2}=-\displaystyle \frac{1-\sqrt{1-{g}^{2}}}{g},{\alpha }_{0,-}=\kappa -\displaystyle \frac{N}{2},\\ & & {E}_{-,N}=-1+N\sqrt{1-{g}^{2}}.\end{array}\end{eqnarray}$The expansion coefficients $a={[{a}_{1},{a}_{2},\cdot \cdot \cdot ,{a}_{N}]}^{T}$ and $b={[{b}_{0},{b}_{1},\cdot \cdot \cdot ,{b}_{N-1}]}^{T}$ are obtained from the recurrence relation$\begin{eqnarray}{R}_{n+1}^{+}{a}_{n+1}+{P}_{n}^{+}{a}_{n}+\displaystyle \frac{{\rm{\Delta }}}{g\xi }{b}_{n}=0,\end{eqnarray}$$\begin{eqnarray}{P}_{n}^{-}{b}_{n}-\displaystyle \frac{{\rm{\Delta }}}{g\xi }{a}_{n}+{Q}_{n-1}^{-}{b}_{n-1}=0,\end{eqnarray}$with ${a}_{0}=0$. For brevity, we present the matrix form$\begin{eqnarray}M\left(\begin{array}{c}a\\ b\end{array}\right)=0,\end{eqnarray}$with$\begin{eqnarray}M=\left(\begin{array}{cc}{M}_{11} & {M}_{12}\\ {M}_{21} & {M}_{22}\end{array}\right).\end{eqnarray}$Here the elements of the matrix M are given as$\begin{eqnarray}\begin{array}{rcl}{M}_{11} & = & \left(\begin{array}{ccccc}{R}_{1}^{+} & & & & \\ {P}_{1}^{+} & {R}_{2}^{+} & & & \\ & {P}_{2}^{+} & \ddots & & \\ & & \ddots & \ddots & \\ & & & {P}_{N-1}^{+} & {R}_{N}^{+}\end{array}\right),\\ {M}_{12} & = & \left(\begin{array}{cccc}\tfrac{{\rm{\Delta }}}{g\xi } & & & \\ & \tfrac{{\rm{\Delta }}}{g\xi } & & \\ & & \ddots & \\ & & & \tfrac{{\rm{\Delta }}}{g\xi }\end{array}\right),\end{array}\end{eqnarray}$$\begin{eqnarray}\begin{array}{rcl}{M}_{21} & = & \left(\begin{array}{cccc}-\tfrac{{\rm{\Delta }}}{g\xi } & & & \\ & -\tfrac{{\rm{\Delta }}}{g\xi } & & \\ & & \ddots & \\ & & & -\tfrac{{\rm{\Delta }}}{g\xi }\end{array}\right),\\ {M}_{22} & = & \left(\begin{array}{ccccc}{Q}_{0}^{-} & {P}_{1}^{-} & & & \\ & {Q}_{1}^{-} & {P}_{2}^{-} & & \\ & & \ddots & \ddots & \\ & & & \ddots & {P}_{N-1}^{-}\\ & & & & {Q}_{N-1}^{-}\end{array}\right).\end{array}\end{eqnarray}$For non-zero values for the expansion coefficients a and b, we further require that the determinant of the matrix M must vanish, thereby leading to the parameter relations for the exact results. For example, in the simplest case of $N=1,2,3$, we get the parameter relations$\begin{eqnarray}N=1,({\rm{\Delta }}-2\kappa +1)({\rm{\Delta }}+2\kappa -1)=0,\end{eqnarray}$$\begin{eqnarray}\begin{array}{rcl}N & = & 2,{{\rm{\Delta }}}^{4}+4({g}^{2}-1+2(1-\kappa )\kappa ){{\rm{\Delta }}}^{2}\\ & & +\,16{\left(1-\kappa \right)}^{2}{\kappa }^{2}=0,\end{array}\end{eqnarray}$$\begin{eqnarray}\begin{array}{l}N=3,{{\rm{\Delta }}}^{3}\pm \left(1-2\kappa \right){{\rm{\Delta }}}^{2}+(8{g}^{2}-4-{\left(1-2\kappa \right)}^{2}){\rm{\Delta }}\\ \quad \pm \,(4{\kappa }^{2}-1)(2\kappa -3)=0.\end{array}\end{eqnarray}$Under the general parameter conditions, the expansion is infinite. One can still choose the following parameters $\eta =-(1-\sqrt{1-{g}^{2}})/g$ and $\xi =-2(1-\sqrt{1-{g}^{2}})/g$, so that ${\lambda }_{3}^{-}=0$ and $1+{\lambda }_{1}^{+}+{\lambda }_{3}^{+}=0$. If we further take the initial conditions a0=0 and b0=1, we require $\alpha =\kappa +(E-1)/2\sqrt{1-{g}^{2}}$. Therefore, we have the three-term recurrence relations for an and bn$\begin{eqnarray}\begin{array}{l}{R}_{n+1}^{+}{P}_{n}^{-}{a}_{n+1}+\left({R}_{n}^{+}{Q}_{n-1}^{-}+{P}_{n}^{+}{P}_{n}^{-}+\displaystyle \frac{{{\rm{\Delta }}}^{2}}{{g}^{2}{\xi }^{2}}\right){a}_{n}\\ \quad +\,{P}_{n-1}^{+}{Q}_{n-1}^{-}{a}_{n-1}=0,\end{array}\end{eqnarray}$$\begin{eqnarray}\begin{array}{l}{R}_{n+1}^{+}{P}_{n+1}^{-}{b}_{n+1}+\left({R}_{n+1}^{+}{Q}_{n}^{-}+{P}_{n}^{+}{P}_{n}^{-}+\displaystyle \frac{{{\rm{\Delta }}}^{2}}{{g}^{2}{\xi }^{2}}\right){b}_{n}\\ \quad +\,{P}_{n}^{+}{Q}_{n-1}^{-}{b}_{n-1}=0.\end{array}\end{eqnarray}$
In the second situation with ${\lambda }_{3}^{+}=0$, we have the parameters and the energies$\begin{eqnarray}\begin{array}{rcl}{\eta }_{+} & = & \displaystyle \frac{{\xi }_{+}}{2}=\displaystyle \frac{1-\sqrt{1-{g}^{2}}}{g},{\alpha }_{0,+}=\kappa -\displaystyle \frac{N}{2},\\ & & {E}_{+,N}=-1+N\sqrt{1-{g}^{2}}.\end{array}\end{eqnarray}$The expansion coefficients $a={[{a}_{0},{a}_{2},\cdot \cdot \cdot ,{a}_{N-1}]}^{T}$ and $b\,={[{b}_{1},{b}_{2},\cdot \cdot \cdot ,{b}_{N}]}^{T}$ satisfy the following equation$\begin{eqnarray}M\left(\begin{array}{c}a\\ b\end{array}\right)=0,\end{eqnarray}$with$\begin{eqnarray}M=\left(\begin{array}{cc}{M}_{11} & {M}_{12}\\ {M}_{21} & {M}_{22}\end{array}\right).\end{eqnarray}$Here the elements of the matrix M are given as$\begin{eqnarray}\begin{array}{rcl}{M}_{11} & = & \left(\begin{array}{cccc}-\tfrac{{\rm{\Delta }}}{g\xi } & & & \\ & -\tfrac{{\rm{\Delta }}}{g\xi } & & \\ & & \ddots & \\ & & & -\tfrac{{\rm{\Delta }}}{g\xi }\end{array}\right),\\ {M}_{12} & = & \left(\begin{array}{ccccc}{R}_{1}^{-} & & & & \\ {P}_{1}^{-} & {R}_{2}^{-} & & & \\ & {P}_{2}^{-} & \ddots & & \\ & & \ddots & \ddots & \\ & & & {P}_{N-1}^{-} & {R}_{N}^{-}\end{array}\right),\end{array}\end{eqnarray}$$\begin{eqnarray}\begin{array}{rcl}{M}_{21} & = & \left(\begin{array}{ccccc}{Q}_{0}^{+} & {P}_{1}^{+} & & & \\ & {Q}_{1}^{+} & {P}_{2}^{+} & & \\ & & \ddots & \ddots & \\ & & & \ddots & {P}_{N-1}^{+}\\ & & & & {Q}_{N-1}^{+}\end{array}\right),\\ {M}_{22} & = & \left(\begin{array}{cccc}\tfrac{{\rm{\Delta }}}{g\xi } & & & \\ & \tfrac{{\rm{\Delta }}}{g\xi } & & \\ & & \ddots & \\ & & & \tfrac{{\rm{\Delta }}}{g\xi }\end{array}\right).\end{array}\end{eqnarray}$In this situation with ${\lambda }_{3}^{+}=0$, the requirement for the non-zero values for the expansion coefficients a and b leads to the same parameter relations with the situation with ${\lambda }_{3}^{-}=0$.
Appendix C. Solutions at the resonant conditions
In this section, we present the solutions at the resonant conditions of ${\rm{\Delta }}=(2\kappa -1)$. Following the results in appendix B, in the case of ${\lambda }_{3}^{-}=0$, we have the exact solutions$\begin{eqnarray}\begin{array}{rcl}{f}_{1}^{-}(z) & = & {a}_{1}^{-}{{\rm{e}}}^{(1-\sqrt{1-{g}^{2}})z/g}\\ & & \times \,F\left(\displaystyle \frac{1}{2}+\kappa ,2\kappa ,-\displaystyle \frac{2-2\sqrt{1-{g}^{2}}}{g}z\right),\\ {f}_{2}^{-}(z) & = & {b}_{0}^{-}{{\rm{e}}}^{(1-\sqrt{1-{g}^{2}})z/g}\\ & & \times \,F\left(-\displaystyle \frac{1}{2}+\kappa ,2\kappa ,-\displaystyle \frac{2-2\sqrt{1-{g}^{2}}}{g}z\right).\end{array}\end{eqnarray}$Here ${a}_{1}^{-}=-{b}_{0}^{-}$. In the case of ${\lambda }_{3}^{+}=0$, we have the exact solutions$\begin{eqnarray}\begin{array}{rcl}{f}_{1}^{+}(z) & = & {a}_{0}^{+}{{\rm{e}}}^{-(1-\sqrt{1-{g}^{2}})z/g}\\ & & \times \,F\left(-\displaystyle \frac{1}{2}+\kappa ,2\kappa ,\displaystyle \frac{2-2\sqrt{1-{g}^{2}}}{g}z\right),\\ {f}_{2}^{+}(z) & = & {b}_{1}^{+}{{\rm{e}}}^{-(1-\sqrt{1-{g}^{2}})z/g}\\ & & \times \,F\left(\displaystyle \frac{1}{2}+\kappa ,2\kappa ,\displaystyle \frac{2-2\sqrt{1-{g}^{2}}}{g}z\right).\end{array}\end{eqnarray}$Here ${b}_{1}^{+}=-{a}_{0}^{+}$. From the relations for the CHF, $F(\alpha ,\gamma ,z)={{\rm{e}}}^{z}F(\gamma -\alpha ,\gamma ,-z)$, it follows that the two eigenfunctions $({f}_{1}^{-}(z),{f}_{2}^{-}(z))$ and $({f}_{1}^{+}(z),{f}_{2}^{+}(z))$ are linearly dependent, and represent the same states. This exact energies are non-degenerate.
Appendix D. Application in the two-photon QRM
In this section, we show that this method can be applied to the two-photon QRM described by the following Hamiltonian (${\hslash }=1$)$\begin{eqnarray}H=\omega {a}^{\dagger }a+{\rm{\Delta }}{\sigma }_{z}+\epsilon {\sigma }_{x}+g{\sigma }_{x}({a}^{\dagger 2}+{a}^{2}).\end{eqnarray}$Here we introduce an additional coupling term $\epsilon {\sigma }_{x}$. In the Fock–Bargmann representation, the bosonic creation and annihilation operators have the form$\begin{eqnarray}a\to \displaystyle \frac{{\rm{d}}}{{\rm{d}}z},\quad {a}^{\dagger }\to z.\end{eqnarray}$The eigenstate $\left|\psi \right\rangle $ can be expressed as$\begin{eqnarray}\left|\psi \right\rangle =\left(\begin{array}{c}{\psi }_{1}(z)\\ {\psi }_{2}(z)\end{array}\right).\end{eqnarray}$From the Schrödinger equation $H| \psi \rangle =E| \psi \rangle $, with the linear combinations, ${\varphi }_{1}=({\psi }_{1}+{\psi }_{2})/2$ and ${\varphi }_{2}=({\psi }_{1}-{\psi }_{2})/2$, we have the two coupled equations for the two component wave functions ${\varphi }_{\mathrm{1,2}}(z)$$\begin{eqnarray}g\displaystyle \frac{{{\rm{d}}}^{2}{\varphi }_{1}}{{\rm{d}}{z}^{2}}+z\displaystyle \frac{{\rm{d}}{\varphi }_{1}}{{\rm{d}}z}+({{gz}}^{2}+\epsilon -E){\varphi }_{1}+{\rm{\Delta }}{\varphi }_{2}=0,\end{eqnarray}$$\begin{eqnarray}g\displaystyle \frac{{{\rm{d}}}^{2}{\varphi }_{2}}{{\rm{d}}{z}^{2}}-z\displaystyle \frac{{\rm{d}}{\varphi }_{2}}{{\rm{d}}z}+({{gz}}^{2}+\epsilon +E){\varphi }_{2}-{\rm{\Delta }}{\varphi }_{1}=0.\end{eqnarray}$Here we have set $\omega =1$. With the following transformations$\begin{eqnarray}{\psi }_{\mathrm{1,2}}(z)={{\rm{e}}}^{-\eta {z}^{2}}{\phi }_{\mathrm{1,2}}(y),y=\xi {z}^{2},\end{eqnarray}$we obtain the similar equations with equations (11) and (12)$\begin{eqnarray}y\displaystyle \frac{{{\rm{d}}}^{2}{\phi }_{1}}{{\rm{d}}{y}^{2}}+\left(\displaystyle \frac{1}{2}+{\lambda }_{1}^{+}y\right)\displaystyle \frac{{\rm{d}}{\phi }_{1}}{{\rm{d}}y}+({\lambda }_{2}^{+}+{\lambda }_{3}^{+}y){\phi }_{1}+\displaystyle \frac{{\rm{\Delta }}}{4g\xi }{\phi }_{2}=0,\end{eqnarray}$$\begin{eqnarray}y\displaystyle \frac{{{\rm{d}}}^{2}{\phi }_{2}}{{\rm{d}}{y}^{2}}+\left(\displaystyle \frac{1}{2}+{\lambda }_{1}^{-}y\right)\displaystyle \frac{{\rm{d}}{\phi }_{2}}{{\rm{d}}y}+({\lambda }_{2}^{-}+{\lambda }_{3}^{-}y){\phi }_{2}-\displaystyle \frac{{\rm{\Delta }}}{4g\xi }{\phi }_{1}=0\end{eqnarray}$with$\begin{eqnarray}{\lambda }_{1}^{\pm }=\displaystyle \frac{-4\eta g\pm 1}{2g\xi },\end{eqnarray}$$\begin{eqnarray}{\lambda }_{2}^{\pm }=\displaystyle \frac{-2\eta g+\epsilon \mp E}{4g\xi },\end{eqnarray}$$\begin{eqnarray}{\lambda }_{3}^{\pm }=\displaystyle \frac{4{\eta }^{2}+g\mp 2\eta }{4g{\xi }^{2}}.\end{eqnarray}$It is evident that our analytical method can be applied directly. Since ${\varphi }_{\mathrm{1,2}}(z)$ have the parity symmetry, ${\varphi }_{\mathrm{1,2}}(-z)\,=\pm {\varphi }_{\mathrm{1,2}}(z)$, one may seek the solutions with the forms of the even-parity and odd-parity functions of z. In the series expansion with the CHFs, the even-parity and odd-parity solutions are given as$\begin{eqnarray}{\phi }_{1}(y)=\sum _{n=0}^{\infty }{a}_{n}\ F({\alpha }_{n},\gamma ,y),\ \ {\phi }_{2}(y)=\sum _{n=0}^{\infty }{b}_{n}F({\alpha }_{n},\gamma ,y),\end{eqnarray}$and$\begin{eqnarray}{\phi }_{1}(y)=\sum _{n=0}^{\infty }{a}_{n}{zF}({\alpha }_{n},\gamma ,y),\ \ {\phi }_{2}(y)=\sum _{n=0}^{\infty }{b}_{n}{zF}({\alpha }_{n},\gamma ,y),\end{eqnarray}$where ${\alpha }_{n}={\alpha }_{0}+n$. By applying our method, for both the even-parity and the odd-parity eigenfunctions, we have ${\eta }_{\pm }={\xi }_{\pm }/2=\pm (\omega -\sqrt{{\omega }^{2}-4{g}^{2}})/4g$. $\gamma =1/2$ and γ=3/2 correspond to the even-parity and odd-parity eigenfunctions, respectively. The different choice of α0 gives the different types of the exact energies. With the choice of ${\alpha }_{0}=-N$, we have the exact energies$\begin{eqnarray}\displaystyle \frac{{E}_{N}^{{\rm{I}},e}}{\omega }=-\displaystyle \frac{1}{2}+\left(2N+\displaystyle \frac{1}{2}\right)\sqrt{1-4\displaystyle \frac{{g}^{2}}{{\omega }^{2}}}\pm \displaystyle \frac{\epsilon }{\omega },\end{eqnarray}$for the even-parity eigenfunctions, and$\begin{eqnarray}\displaystyle \frac{{E}_{N}^{{\rm{I}},o}}{\omega }=-\displaystyle \frac{1}{2}+\left(2N+\displaystyle \frac{3}{2}\right)\sqrt{1-4\displaystyle \frac{{g}^{2}}{{\omega }^{2}}}\pm \displaystyle \frac{\epsilon }{\omega },\end{eqnarray}$for the odd-parity eigenfunctions. These exact results indicate that due to the presence of the additional term εσx, these exact energies associated with polynomial eigenfunctions do not correspond to the level-crossings. If we take ${\alpha }_{0}=(1-2N\pm 2\epsilon /\sqrt{{\omega }^{2}-4{g}^{2}})/4$ and ${\alpha }_{0}=(3-2N\pm 2\epsilon /\sqrt{{\omega }^{2}-4{g}^{2}})/4$ for the even-parity and the odd-parity eigenfunctions, we have the exact energies$\begin{eqnarray}\displaystyle \frac{{E}_{N}^{\mathrm{II},e,o}}{\omega }=-\displaystyle \frac{1}{2}+N\sqrt{1-4\displaystyle \frac{{g}^{2}}{{\omega }^{2}}}.\end{eqnarray}$Here we have N≥1. In the simplest case of N=1, the parameter relation is given as$\begin{eqnarray}\displaystyle \frac{{{\rm{\Delta }}}^{2}}{{\omega }^{2}}-\displaystyle \frac{1}{4}+\displaystyle \frac{4}{1-4{g}^{2}/{\omega }^{2}}\displaystyle \frac{{\epsilon }^{2}}{{\omega }^{2}}=0.\end{eqnarray}$This means that for the resonance condition of 2Δ=ω with ε=0, we have the exact energy ${E}_{N}^{\mathrm{II},e,o}/\omega \,=-1/2+\sqrt{1-4{g}^{2}/{\omega }^{2}}$.