# Analytical approximations for generalized quantum Rabi models

Chon-Fai Kam  Email<EMAIL_ADDRESS>Department of Physics, University
at Buffalo, SUNY, Buffalo, New York 14260, USA  Yang Chen  Email:
<EMAIL_ADDRESS>Department of Mathematics, Faculty of Science and
Technology, University of Macau, Avenida da Universidade, Taipa, Macau, China

###### Abstract

The quantum Rabi model is essential for understanding interacting quantum
systems. It serves as the simplest non-integrable yet solvable model
describing the interaction between a two-level system and a single mode of a
bosonic field. In this study, we delve into the exploration of the generalized
quantum Rabi model, wherein the bosonic mode of the field undergoes squeezing.
Utilizing the Segal-Bargmann representation of the infinite-dimensional
Hilbert space, we demonstrate that the energy spectrum of the generalized
quantum Rabi model, when both the Rabi coupling strength and the squeezing
strength are not significantly large compared to the field mode frequency, can
be analytically determined by a bi-confluent Fuchsian equation with two
regular singularities at $0$ and $1$ and an irregular singularity of rank two
at infinity.

##  I Introduction

The quantum Rabi model is probably one of the simplest physical systems
surpassing the quantum harmonic oscillator. It describes the interaction of a
single two-level system and a single bosonic mode. Rabi’s original paper
concerned the effect of a rapidly varying magnetic field on an oriented atom
possessing nuclear spin rabi1936process ; rabi1937space . Over the last eighty
years, it has become the cornerstone for theoretically explaining various
physical systems, including cavity quantum electrodynamics (QED)
gleyzes2007quantum ; guerlin2007progressive , superconducting circuits
forn2010observation , mechanical oscillators lahaye2009nanomechanical , and
photonic crystals liu2017quantum , among others. Although the quantum Rabi
model possesses a degree of freedom of two, which is smaller than that of a
hydrogen atom, the energy spectrum of the infinite-dimensional state space has
only been obtained through numerical diagonalization of the truncated finite-
dimensional Hilbert space for an extended period. It was not until 2011 that
Braak demonstrated that the exact solvability of the quantum Rabi model
braak2011integrability . Braak’s approach is based on the Segal-Bargmann
holomorphic representation of Hilbert space bargmann1961hilbert ;
segal1963mathematical . In this representation, the state vector is expressed
as an integral transform known as the Segal-Bargmann transform. The kernel of
this transform is the inner product between a Gilmore-Perelomov coherent state
zhang1990coherent ; kam2023coherent  and the state vector itself. In Braak’s
original paper, the regular part of the spectrum is determined by the zeros of
a transcendental $G$-function associated with the Hamiltonian. This
$G$-function serves as a generalized spectral determinant with only zeros on
the real axis, comprising all regular subsets of the discrete eigenvalues.

Following Braak’s groundbreaking work, it was further recognized that the
holomorphic function in Segal-Bargmann space, satisfying the square-
integrability condition, is governed by the confluent Heun equation
zhong2013analytical . This equation, featuring two regular singularities at
$0$ and $1$ and an irregular singularity at infinity, represents the confluent
form of the general Heun equation ronveaux1995heun . The latter includes four
regular singularities in the extended complex plane, representing the simplest
differential equation that cannot be solved by conventional special functions.
The quantization condition is determined by gluing the local holomorphic Heun
functions near each regular singularity into a global entire function on the
complex plane. The regular spectrum of the quantum Rabi model is determined by
a transcendental entire function constructed by the Wronskian of the local
Heun functions. The exceptional eigenvalues in the energy spectrum, called
Judd points judd1979exact , are determined by the condition that the confluent
Heun function is truncated into the confluent Heun polynomial, which is
regular at all singularities on whole complex plane. In this regard, it is no
wonder that the exact solution of the quantum Rabi model has not been
discovered for many decades. All previously solvable models, such as the
quantum harmonic oscillator, the hydrogen atom, and the symmetric top,
incorporate the hypergeometric function and its confluent form pauli1973pauli
. These functions are simpler than Heun functions in the sense that they have
at most three regular singularities in the extended complex plane. Besides,
the coefficients of the local Heun series are determined by a three-term
recursion relation, in contrast to the conventional two-term recursion
relation found in all hypergeometric series ronveaux1995heun . The versatility
of Heun functions in solving complicated differential equations makes them a
valuable tool in theoretical physics. In recent years, Heun functions and
their confluent forms have found numerous applications, including solving the
time-dependent quantum two-state problem ishkhanyan2014fifteen ;
ishkhanyan2015thirty ; kam2020analytical ; kam2021analytical , describing the
quantum Ising chain in a periodic transverse magnetic field dorosz2008work ,
analyzing two-dimensional massless Dirac fermions in spatially inhomogeneous
fields dorosz2008work , and addressing non-Hermitian physics in optical
waveguides kam2021non ; kam2023analytical , among others.

In this study, we investigate the energy spectrum of a generalized quantum
Rabi model, where the bosonic mode experiences squeezing. This model can be
employed to describe the interaction of a two-level atom with squeezed light
gerry1988two , or characterize a quantum dot in a squeezed cavity field
del2010two . Here, squeezed states refer to a class of nonclassical states
characterized by modified quantum noise profiles zhang1990coherent . By
introducing quantum correlations among photons, quantum fluctuations are
periodically reduced to below the standard quantum limit in one field
quadrature while simultaneously increased in the conjugate quadrature
kam2023coherent . Squeezed states have been thoroughly investigated across
diverse research domains. Notably, high-intensity squeezed light has found
application in cutting-edge gravitational wave detectors tse2019quantum ;
acernese2019increasing , leading to a nearly tenfold improvement in
sensitivity dwyer2022squeezing . Within the realm of quantum information
processing, squeezed states have been utilized in continuous-variable quantum
key distribution gehring2015implementation , quantum sensing lawrie2019quantum
, and high-precision cavity spectroscopy junker2021high . Additionally, in the
field of semiconductor quantum dots, squeezed states have been proposed to
enhance the signal-to-noise ratio (SNR) and improve the readout fidelity of
the output signal within sub-microsecond timescales kam2023sub ; kam2024fast .

For the generalized quantum Rabi model with squeezing, we derive the fourth-
order differential equation that governs the holomorphic function in the
Segal-Bargmann space. Subsequently, we obtain its series on the complex plane,
whose coefficients are determined by a nine-term recursion relation. When both
the Rabi coupling strength and the squeezing strength are not large compared
to the field mode frequency, the fourth-order differential equation that
governs the holomorphic function in the Segal-Bargmann space can be truncated
to a second-order one. The approximated equation is the bi-confluent form of a
general Fuchsian equation with five singularities in the extended complex
plane, whose coefficients are determined by a four-term recursion relation. As
a result, the regular parts of the spectrum are determined by a transcendental
$G$-function, which is the Wronskian of the local holomorphic series solutions
at the two regular singularities $0$ and $1$. Additionally, the exceptional
eigenvalues in the energy spectrum, namely, the Judd points, are determined by
the condition that the bi-confluent Fuchsian function is truncated into a
polynomial.

##  II The generalized quantum Rabi model

Let us consider the generalized Rabi model where the Hamiltonian is described
by

$H=\Delta\sigma_{z}+\epsilon\sigma_{x}+\omega
a^{\dagger}a+g\sigma_{x}(a^{\dagger}+a)+\lambda\sigma_{x}(a^{\dagger
2}+a^{2}),$ (1)

where $a$ and $a^{\dagger}$ are the annihilation and creation operators for a
single bosonic mode with frequency $\omega$, $\sigma_{x}$ and $\sigma_{z}$ are
the Pauli matrices for a two-level system with level splitting 2$\Delta$, the
term $\epsilon\sigma_{x}$ describes a spontaneous transition of the two-level
system which is not driven by the field, the term proportional to $g$
signifies a direct interaction between the two systems, and the term
proportional to $\lambda$ introduces squeezing into the bosonic mode. The
Hamiltonian (1) is a single-mode spin boson model. It reduces to the standard
quantum Rabi model when $\epsilon=0$ and $\lambda=0$. In contrast, the
Hamiltonian (1) reduces to the two-photon Rabi model when $\epsilon=0$ and
$g=0$, and it reduces to the asymmetric quantum Rabi model when $\lambda=0$.
To solve the generalized Rabi model, one can employ the continuous
representation of the two-photon group in the Segal-Bargmann space of
holomorphic functions in a complex variable $z$. In this space, the generators
of the two-photon group have the form

$a^{\dagger 2}\rightarrow
z^{2},a^{2}\rightarrow\frac{d^{2}}{dz^{2}},a^{\dagger}a\rightarrow
z\frac{d}{dz},a\rightarrow\frac{d}{dz},a^{\dagger}\rightarrow z.$ (2)

To study the energy spectrum of the generalized Rabi model, one can write the
eigenstate $|\psi\rangle$ of the Hamiltonian (1) as a two-component wave
function $(\psi_{1},\psi_{2})^{\top}$. From the Schrödinger equation
$H|\psi\rangle=E|\psi\rangle$, the Hamiltonian becomes

$H=\begin{pmatrix}\omega
a^{\dagger}a+\Delta&g(a+a^{\dagger})+\lambda(a^{2}+a^{\dagger 2})+\epsilon\\\
g(a+a^{\dagger})+\lambda(a^{2}+a^{\dagger 2})+\epsilon&\omega
a^{\dagger}a-\Delta\end{pmatrix}.$ (3)

To better study the eigenvalue problem, one can introduce a new two-component
wave function $(\phi_{1},\phi_{2})^{\top}$ via

$\begin{pmatrix}\psi_{1}\\\ \psi_{2}\end{pmatrix}=U\begin{pmatrix}\phi_{1}\\\
\phi_{2}\end{pmatrix},U\equiv\frac{1}{\sqrt{2}}(\sigma_{z}+\sigma_{x}).$ (4)

Using the properties $U^{-1}\sigma_{z}U=\sigma_{x}$ and
$U^{-1}\sigma_{x}U=\sigma_{z}$, the Hamiltonian $H^{\prime}\equiv U^{-1}HU$ in
the new basis becomes

$H^{\prime}=\epsilon\sigma_{z}+\Delta\sigma_{x}+\omega
a^{\dagger}a+g\sigma_{z}(a^{\dagger}+a)+\lambda\sigma_{z}(a^{\dagger
2}+a^{2}).$ (5)

In the Segal-Bargmann space of holomorphic functions, the generalized Rabi
model can be mapped into a set of two coupled second-order ordinary
differential equations for the wave function components $\phi_{1}$ and
$\phi_{2}$

$\displaystyle\lambda\frac{d^{2}\phi_{1}}{dz^{2}}+(g+\omega
z)\frac{d\phi_{1}}{dz}+(\epsilon-E+gz+\lambda z^{2})\phi_{1}$
$\displaystyle=-\Delta\phi_{2},$ (6a)
$\displaystyle\lambda\frac{d^{2}\phi_{2}}{dz^{2}}+(g-\omega
z)\frac{d\phi_{2}}{dz}+(\epsilon+E+gz+\lambda z^{2})\phi_{2}$
$\displaystyle=\Delta\phi_{1}.$ (6b)

Eq. (6b) becomes Eq. (6a) under the transformation
$\epsilon\rightarrow-\epsilon$, $g\rightarrow-g$ and
$\lambda\rightarrow-\lambda$. To solve the coupled second-order differential
equation, let us introduce the following differential operators
$D_{z}\equiv\lambda\frac{d^{2}}{dz^{2}}+c_{1}(z)\frac{d}{dz}+c_{2}(z)$, and
$\bar{D}_{z}\equiv\lambda\frac{d^{2}}{dz^{2}}+\bar{c}_{1}(z)\frac{d}{dz}+\bar{c}_{2}(z)$,
where $c_{1}(z)\equiv g+\omega z$, $c_{2}(z)\equiv\epsilon-E+gz+\lambda
z^{2}$, $\bar{c}_{1}(z)\equiv g-\omega z$, and
$\bar{c}_{2}(z)\equiv\epsilon+E+gz+\lambda z^{2}$. Then the wave amplitude
$\phi_{1}$ is governed by $(\bar{D}_{z}D_{z}+\Delta^{2})\phi_{1}$, or
equivalently

$\displaystyle\lambda^{2}\phi_{1}^{(4)}+\lambda(c_{1}+\bar{c}_{1})\phi_{1}^{(3)}+(2\lambda
c_{1}^{\prime}+c_{1}\bar{c}_{1}+\lambda\bar{c}_{2})\phi_{1}^{\prime\prime}$
$\displaystyle+$ $\displaystyle(\lambda
c_{1}^{\prime\prime}+\bar{c}_{1}c_{1}^{\prime}+c_{1}\bar{c}_{2})\phi_{1}^{\prime}+(\lambda
c_{2}^{\prime\prime}+\bar{c}_{1}c_{2}^{\prime}+\bar{c}_{2}c_{2}+\Delta^{2})\phi_{1}=0,$

where $c_{j}^{\prime}$ denotes the derivative of $c_{j}$ with respect to $z$.
A direct computation yields the following fourth-order differential equation

$\displaystyle\lambda^{2}\phi_{1}^{(4)}+2\lambda
g\phi_{1}^{(3)}+[g^{2}-\omega^{2}z^{2}+\lambda(2\omega+\epsilon+E+gz+\lambda
z^{2})]\phi_{1}^{\prime\prime}$ $\displaystyle+$
$\displaystyle[\omega(g-\omega z)+(g+\omega z)(\epsilon+E+gz+\lambda
z^{2})]\phi_{1}^{\prime}+[2\lambda^{2}$ $\displaystyle+$
$\displaystyle(g-\omega z)(g+2\lambda z)+(\epsilon+gz+\lambda
z^{2})^{2}-E^{2}+\Delta^{2}]\phi_{1}=0.$ (7)

###  II.1 The uncoupled case with $\Delta=0$

For the special case when $\Delta=0$, $\phi_{1}$ and $\phi_{2}$ are governed
by two independent second-order differential equations

$\displaystyle\lambda\frac{d^{2}\phi_{1}}{dz^{2}}+(g+\omega
z)\frac{d\phi_{1}}{dz}+(\epsilon-E+gz+\lambda z^{2})\phi_{1}$
$\displaystyle=0,$ (8a)
$\displaystyle\lambda\frac{d^{2}\phi_{2}}{dz^{2}}+(g-\omega
z)\frac{d\phi_{2}}{dz}+(\epsilon+E+gz+\lambda z^{2})\phi_{2}$
$\displaystyle=0,$ (8b)

For a general second-order ordinary differential equation

$\frac{d^{2}\phi}{dz^{2}}+p(z)\frac{d\phi}{dz}+q(z)\phi=0,$ (9)

one can introduce the following transformation

$\phi(z)=u(z)\exp\left(-\frac{1}{2}\int_{0}^{z}p(z)dz^{\prime}\right),$ (10)

to obtain its second canonical form

$\displaystyle\frac{d^{2}u}{dz^{2}}+Q(z)u=0,Q(z)\equiv-\frac{1}{2}p^{\prime}(z)-\frac{1}{4}p^{2}(z)+q(z).$
(11)

Hence, using the transformations
$\phi_{1}(z)=u_{1}(z)\exp(-\frac{g}{2\lambda}z-\frac{\omega}{4\lambda}z^{2})$,
one obtains the Weber equation

$\frac{d^{2}u_{1}}{d\zeta_{1}^{2}}=\left(\frac{1}{4}\zeta_{1}^{2}+a_{1}\right)u_{1},$
(12)

where

$\displaystyle\zeta_{1}$
$\displaystyle\equiv\left(\frac{\omega^{2}}{\lambda^{2}}-4\right)^{\frac{1}{4}}\left(z+\frac{g}{\omega+2\lambda}\right),$
(13a) $\displaystyle a_{1}$
$\displaystyle\equiv\frac{1}{\lambda}\left(\frac{\omega^{2}}{\lambda^{2}}-4\right)^{-\frac{1}{2}}\left(E+\frac{g^{2}}{\omega+2\lambda}+\frac{\omega}{2}-\epsilon\right).$
(13b)

The function $u_{1}(\zeta_{1})$ in Eqs. (12) can be solved by the parabolic
cylinder functions. Explicitly, the two independent even and odd solutions of
the function $u_{1}(\zeta_{1})$ are denoted by $U_{e}(a_{1},\zeta_{1})$ and
$U_{o}(a_{1},\zeta_{1})$ via

$\displaystyle U_{e}(a_{1},\zeta_{1})$
$\displaystyle\equiv\>\>\>\>\>\exp\left(-\frac{\zeta_{1}^{2}}{4}\right){}_{1}F_{1}\left(\frac{a_{1}}{2}+\frac{1}{4};\frac{1}{2};\frac{\zeta_{1}^{2}}{2}\right),$
(14a) $\displaystyle U_{o}(a_{1},\zeta_{1})$
$\displaystyle\equiv\zeta_{1}\exp\left(-\frac{\zeta_{1}^{2}}{4}\right){}_{1}F_{1}\left(\frac{a_{1}}{2}+\frac{3}{4};\frac{3}{2};\frac{\zeta_{1}^{2}}{2}\right),$
(14b)

where ${}_{1}F_{1}(a;b;\zeta_{1})$ is the Kummer confluent hypergeometric
function. From the Weber equation Eq. (12), the energy spectrum of the
generalized Rabi model with $\Delta=0$ has two branches, where the positive
branch with $\sigma_{z}=1$ is determined by the condition
$a_{1}=n_{1}+\frac{1}{2}$, or equivalently

$E_{n_{1}}=\sqrt{\omega^{2}-4\lambda^{2}}\left(n_{1}+\frac{1}{2}\right)-\frac{g^{2}}{\omega+2\lambda}-\frac{\omega}{2}+\epsilon,$
(15)

and the negative branch of the energy spectrum with $\sigma_{z}=-1$ is
determined by using the transformation $\epsilon\rightarrow-\epsilon$,
$g\rightarrow-g$ and $\lambda\rightarrow-\lambda$, which yields

$E_{n_{2}}=\sqrt{\omega^{2}-4\lambda^{2}}\left(n_{2}+\frac{1}{2}\right)-\frac{g^{2}}{\omega-2\lambda}-\frac{\omega}{2}-\epsilon.$
(16)

###  II.2 The asymmetric Rabi model with $\lambda=0$

For the special case with $\lambda=0$, which is associated with the asymmetric
Rabi model, the fourth-order differential equation Eq. (II) is simplified to a
second-order differential equation

$\displaystyle(g-\omega z)(g+\omega
z)\phi_{1}^{\prime\prime}+\left[\omega(g-\omega z)+(g+\omega
z)(\epsilon+E+gz)\right]\phi_{1}^{\prime}$ $\displaystyle+[(g-\omega
z)g+(\epsilon+gz)^{2}-E^{2}+\Delta^{2}]\phi_{1}=0.$ (17)

A direct computation yields

$\phi_{1}^{\prime\prime}+\left(A_{1}+\frac{A_{2}}{z-q}+\frac{A_{3}}{z+q}\right)\phi_{1}^{\prime}+\left(B_{1}+\frac{B_{2}}{z-q}+\frac{B_{3}}{z+q}\right)\phi_{1}=0,$
(18)

where $q\equiv g/\omega$ and

$\displaystyle A_{1}$
$\displaystyle\equiv-q,B_{1}\equiv-q^{2},A_{2}\equiv-q^{2}-\frac{\epsilon+E}{\omega},A_{3}\equiv
1,$ $\displaystyle B_{2}$ $\displaystyle\equiv-\frac{q^{3}}{2}-\frac{\epsilon
q}{\omega}-\frac{\epsilon^{2}-E^{2}+\Delta^{2}}{2\omega g},$ $\displaystyle
B_{3}$ $\displaystyle\equiv\>\>\>\frac{q^{3}}{2}-\frac{\epsilon
q}{\omega}+\frac{\epsilon^{2}-E^{2}+\Delta^{2}}{2\omega g}+q.$ (19)

After the linear variable transformation $\zeta\equiv(\omega z+g)/2g$, one
obtains

$\frac{d^{2}\phi_{1}}{d\zeta^{2}}+\left(\alpha_{1}+\frac{\alpha_{2}}{\zeta-1}+\frac{\alpha_{3}}{\zeta}\right)\frac{d\phi_{1}}{d\zeta}+\left(\beta_{1}+\frac{\beta_{2}}{\zeta-1}+\frac{\beta_{3}}{\zeta}\right)\phi_{1}=0,$
(20)

where

$\displaystyle\alpha_{1}$
$\displaystyle\equiv-2q^{2},\alpha_{2}=-q^{2}-\frac{\epsilon+E}{\omega},\alpha_{3}=1,\beta_{1}\equiv-4q^{4},$
$\displaystyle\beta_{2}$ $\displaystyle\equiv-q^{4}-\frac{2\epsilon
q^{2}}{\omega}-\frac{\epsilon^{2}-E^{2}+\Delta^{2}}{\omega^{2}},$
$\displaystyle\beta_{3}$ $\displaystyle\equiv\>\>\>q^{4}-\frac{2\epsilon
q^{2}}{\omega}+\frac{\epsilon^{2}-E^{2}+\Delta^{2}}{\omega^{2}}+2q^{2}.$ (21)

Finally, the above equation can be reduced to the standard confluent Heun
equation after the transformation $\phi_{1}(\zeta)\equiv
e^{k\zeta}w_{1}(\zeta)$

$\frac{d^{2}w_{1}}{d\zeta^{2}}+\left(\alpha+\frac{\beta+1}{\zeta}+\frac{\gamma}{\zeta-1}\right)\frac{dw_{1}}{d\zeta}+\left(\frac{\mu}{\zeta}+\frac{\nu}{\zeta-1}\right)w_{1}=0,$
(22)

where $\alpha\equiv\alpha_{1}+2k$, $\beta=0$, $\gamma\equiv\alpha_{2}$,
$\mu\equiv k\alpha_{3}+\beta_{3}$, $\nu\equiv k\alpha_{2}+\beta_{2}$, and $k$
is determined by the relation $k^{2}+\alpha_{1}k+\beta_{1}=0$, which is
explicitly solved by $k_{\pm}=(1\pm\sqrt{5})g^{2}/\omega^{2}$. Eq. (22) is
solved by the confluent Heun function, which has the power series near the
singularity $\zeta=0$

$w_{1}^{(0)}(\zeta)\equiv\sum_{n=0}^{\infty}a_{n}(\alpha,\beta,\gamma,\mu,\nu)\zeta^{n},$
(23)

where the coefficients $a_{n}$ are determined by the three-term recurrence
relation

$a_{n+1}=\frac{n^{2}+(\beta+\gamma-\alpha)-\mu}{(n+1)(n+\beta+1)}a_{n}+\frac{\alpha(n-1)+(\mu+\nu)}{(n+1)(n+\beta+1)}a_{n-1},$
(24)

with initial conditions $a_{-1}=0$ and $a_{0}=1$. This series converges for
$|\zeta|<1$. One can also expand the confluent Heun function near the
singularity $\zeta=1$ into power series

$w_{1}^{(1)}(\zeta)=\sum_{n=0}^{\infty}b_{n}(\alpha,\beta,\gamma,\mu,\nu)(\zeta-1)^{n},$
(25)

where the coefficients $b_{n}$ are determined by the three-term recurrence
relation

$b_{n+1}+\frac{n(\alpha+\beta+\gamma+n)+\nu}{(n+1)(n+\gamma)}b_{n}+\frac{(n-1)\alpha+\mu+\nu}{(n+1)(n+\gamma)}b_{n-1}=0,$
(26)

with initial conditions $b_{-1}=0$ and $b_{0}=1$. This series defines a local
holomorphic function on $|\zeta-1|<1$. To determine the quantization
condition, one needs to analyze the analytic continuation of the local
holomorphic series solutions. The local holomorphic functions
$w_{1}^{(0)}(\zeta)$ and $w_{1}^{(1)}(\zeta)$ at the regular singularities $0$
and $1$ can be extended to a global entire function on the complex plane when
$w_{1}^{(0)}(\zeta)$ and $w_{1}^{(1)}(\zeta)$ and their derivatives satisfy
the condition

$G(E)\equiv
w_{1}^{(0)}(\zeta;E)w_{1}^{(1)\prime}(\zeta;E)-w_{1}^{(1)}(\zeta;E)w_{1}^{(0)\prime}(\zeta;E),$
(27)

where $G(E)$ is the Wronskian of the local holomorphic series at the regular
singularities $0$ and $1$, regarded as a function of $E$, and vanishes on the
quantized energies.

###  II.3 Two-photon quantum Rabi model with $g=0$

We now consider another case when $\lambda\neq 0$ but with $g=0$, which is
associate with the two-photon quantum Rabi model. The wave amplitude
$\phi_{1}$ is thus governed by

$\phi_{1}^{(4)}+(A_{1}+A_{2}z^{2})\phi_{1}^{\prime\prime}+(B_{1}z+B_{2}z^{3})\phi_{1}^{\prime}+(C_{1}+C_{2}z^{2}+z^{4})\phi_{1}=0,$
(28)

where

$\displaystyle A_{1}\equiv\frac{2\omega+\epsilon+E}{\lambda},A_{2}\equiv
1-\frac{\omega^{2}}{\lambda^{2}},B_{1}\equiv\frac{\omega(\epsilon+E)-\omega^{2}}{\lambda^{2}},$
$\displaystyle
B_{2}\equiv\frac{\omega}{\lambda},C_{1}\equiv\frac{\epsilon^{2}-E^{2}-\Delta^{2}}{\lambda^{2}}+2,C_{2}\equiv\frac{2(\epsilon-\omega)}{\lambda^{2}}.$
(29)

Eq. (II.3) is solved by the following series

$\phi_{1}(z)\equiv\sum_{n=0}^{\infty}a_{n}(A_{1},A_{2},B_{1},B_{2},C_{1},C_{2})z^{n},$
(30)

where the coefficients $a_{n}$ are determined by the five-term recurrence
relation

$\displaystyle(n+4)(n+3)(n+2)(n+1)a_{n+4}+(n+2)[(n+1)A_{1}+B_{2}]a_{n+2}$
$\displaystyle+[n(n-1)A_{2}+nB_{1}+C_{1}]a_{n}+C_{2}a_{n-2}+a_{n-4}=0.$ (31)

The above series converges for all finite $|z|$, and thus it defines an entire
function on the complex plane.

###  II.4 The general case with $g\neq 0$ and $\lambda\neq 0$

We now return to the general case where both $g$ and $\lambda$ are non-zero,
where the wave amplitude $\phi_{1}$ is governed by

$\displaystyle\frac{d^{4}\phi_{1}}{dz^{4}}+A_{1}\frac{d^{3}\phi_{1}}{dz^{3}}+(B_{1}+B_{2}z+B_{3}z^{2})\frac{d^{2}\phi_{1}}{dz^{2}}$
$\displaystyle+(C_{1}+C_{2}z+C_{3}z^{2}+C_{4}z^{3})\frac{d\phi_{1}}{dz}$
$\displaystyle+(D_{1}+D_{2}z+D_{3}z^{2}+D_{4}z^{3}+z^{4})\phi_{1}=0,$ (32)

where

$\displaystyle
A_{1}\equiv\frac{2g}{\lambda},B_{1}\equiv\frac{g^{2}}{\lambda^{2}}+\frac{2\omega+\epsilon+E}{\lambda},$
$\displaystyle B_{2}\equiv\frac{g}{\lambda},B_{3}\equiv
1-\frac{\omega^{2}}{\lambda^{2}},C_{1}\equiv\frac{g(\omega+\epsilon+E)}{\lambda^{2}},$
$\displaystyle
C_{2}\equiv\frac{-\omega^{2}+\omega(\epsilon+E)+g^{2}}{\lambda^{2}},C_{3}\equiv\frac{\omega
g}{\lambda^{2}}+\frac{g}{\lambda},$ $\displaystyle
C_{4}\equiv\frac{\omega}{\lambda},D_{1}\equiv
2+\frac{g^{2}+\epsilon^{2}-E^{2}-\Delta^{2}}{\lambda^{2}},$ $\displaystyle
D_{2}\equiv\frac{g(2\epsilon-\omega)}{\lambda^{2}}+\frac{2g}{\lambda},D_{3}\equiv\frac{g^{2}}{\lambda^{2}}+\frac{2(\epsilon-\omega)}{\lambda},D_{4}\equiv\frac{2g}{\lambda}.$
(33)

Eq. (II.4) is solved by the following series

$\phi_{1}(z)\equiv\sum_{n=0}^{\infty}a_{n}(A_{1},B_{1},B_{2},B_{3},C_{1},C_{2},C_{3},C_{4},D_{1},D_{2},D_{3},D_{4})z^{n},$
(34)

where the coefficients $a_{n}$ satisfy the nine-term recurrence relation

$\displaystyle(n+4)(n+3)(n+2)(n+1)a_{n+4}+(n+3)(n+2)(n+1)A_{1}a_{n+3}$
$\displaystyle+$
$\displaystyle(n+2)(n+1)B_{1}a_{n+2}+(n+1)(nB_{2}+C_{1})a_{n+2}$
$\displaystyle+$
$\displaystyle[n(n-1)B_{3}+nC_{2}+D_{1}]a_{n}+[(n-1)C_{3}+D_{2}]a_{n-1}$
$\displaystyle+$
$\displaystyle[(n-2)C_{4}+D_{3}]a_{n-2}+D_{4}a_{n-3}+a_{n-4}=0.$ (35)

Interestingly, when both $g$ and $\lambda$ are small but non zero parameters,
one can neglect all the terms proportional to $\lambda^{2}$ and $\lambda g$
and obtain

$(z^{2}-q^{2})\phi_{1}^{\prime\prime}=(A_{1}+A_{2}z+A_{3}z^{2}+A_{4}z^{3})\phi_{1}^{\prime}+(B_{1}+B_{2}z+B_{3}z^{2})\phi_{1},$
(36)

where

$\displaystyle
q\equiv\sqrt{\frac{g^{2}+\lambda(\epsilon+E)}{\omega^{2}}+\frac{2\lambda}{\omega}},A_{1}\equiv\frac{g(\epsilon+E)}{\omega^{2}}+\frac{g}{\omega},$
$\displaystyle
A_{2}\equiv\frac{g^{2}}{\omega^{2}}+\frac{\epsilon+E}{\omega}-1,A_{3}\equiv\frac{g}{\omega},$
$\displaystyle
A_{4}\equiv\frac{\lambda}{\omega},B_{1}\equiv\frac{g^{2}+\epsilon^{2}-E^{2}-\Delta^{2}}{\omega^{2}},$
$\displaystyle B_{2}\equiv\frac{2\epsilon
g}{\omega^{2}}-\frac{g}{\omega},B_{3}\equiv\frac{g^{2}+2\epsilon\lambda}{\omega^{2}}-\frac{2\lambda}{\omega}.$
(37)

A direct computation yields

$\displaystyle-$
$\displaystyle\frac{d^{2}\phi_{1}}{dz^{2}}+\left(A_{4}z+A_{3}+\frac{\beta_{+}}{z-q}+\frac{\beta_{-}}{z+q}\right)\frac{d\phi_{1}}{dz}$
$\displaystyle+$
$\displaystyle\left(B_{3}+\frac{\gamma_{+}}{z-q}+\frac{\gamma_{-}}{z+q}\right)\phi_{1}=0,$
(38)

where

$\beta_{\pm}\equiv\frac{1}{2}\left[(A_{4}\pm A_{3})q^{2}+A_{2}\pm
A_{1}\right],\gamma_{\pm}\equiv\frac{1}{2}\left[B_{2}\pm\frac{B_{3}q^{2}+B_{1}}{q}\right].$
(39)

After the transformation $\zeta\equiv z/2q+1/2$, Eq. (II.4) becomes

$\frac{d^{2}\phi_{1}}{d\zeta^{2}}=\left(\alpha_{1}\zeta+\alpha_{2}+\frac{\beta_{1}}{\zeta-1}+\frac{\beta_{2}}{\zeta}\right)\frac{d\phi_{1}}{d\zeta}+\left(\gamma_{1}+\frac{\gamma_{2}}{\zeta-1}+\frac{\gamma_{3}}{\zeta}\right)\phi_{1},$
(40)

where $\alpha_{1}\equiv 2qA_{4}$, $\alpha_{2}\equiv 2qA_{3}-2q^{2}A_{4}$,
$\beta_{1}\equiv\beta_{+}$, $\beta_{2}\equiv\beta_{-}$, $\gamma_{1}\equiv
4q^{2}B_{3}$, $\gamma_{2}\equiv 2q\gamma_{+}$ and $\gamma_{3}\equiv
2q\gamma_{-}$. Eq. (40) has three singularities, two regular singularities at
$0$ and $1$, and an irregular singularity of rank $2$ at $\infty$. It is the
Biconfluent form of an Fuchsian equation with five regularities. The local
holomorphic solution around the regular singularity $0$ is determined by the
power series

$\phi_{1}^{(0)}(\zeta)\equiv\sum_{n=0}^{\infty}a_{n}(\alpha_{1},\alpha_{2},\beta_{1},\beta_{2},\gamma_{1},\gamma_{2},\gamma_{3})\zeta^{n}.$
(41)

The coefficients $a_{n}$ satisfy the four-term recurrence relation

$\displaystyle(n+1)(n+\beta_{2})a_{n+1}=[n(n-1)+n\mu-\gamma_{3}]a_{n}$
$\displaystyle+[(n-1)\nu+\gamma]a_{n-1}+[(n-2)\alpha_{1}+\gamma_{1}]a_{n-2},$
(42)

where $\mu\equiv\beta_{1}+\beta_{2}-\alpha_{2}$,
$\nu\equiv\alpha_{2}-\alpha_{1}$,
$\gamma\equiv\gamma_{2}+\gamma_{3}-\gamma_{1}$. Eq. (II.4) reduces to the
three-term recurrence relation for the confluent Heun function if
$\alpha_{1}=\gamma_{1}=0$. It defines a local holomorphic function on
$|\zeta|<1$. Similarly, the local holomorphic solution around the regular
singularity $1$ is determined by the power series

$\phi_{1}^{(1)}(\zeta)\equiv\sum_{n=0}^{\infty}b_{n}(\alpha_{1},\alpha_{2},\beta_{1},\beta_{2},\gamma_{1},\gamma_{2},\gamma_{3})(\zeta-1)^{n}.$
(43)

The coefficients $b_{n}$ satisfy the four-term recurrence relation

$\displaystyle(n+1)(n-\beta_{1})b_{n+1}=[n\delta+\gamma_{2}-n(n-1)]b_{n}$
$\displaystyle+[(n-1)\eta+\kappa]b_{n-1}+[(n-2)\alpha_{2}+\gamma_{1}]b_{n-2},$
(44)

where $\delta\equiv\alpha_{1}+\alpha_{2}+\beta_{1}+\beta_{2}$, $\eta\equiv
2\alpha_{1}+\alpha_{2}$, $\kappa\equiv\gamma_{1}+\gamma_{2}+\gamma_{3}$. It
defines a local holomorphic function on $|\zeta-1|<1$. The local holomorphic
series $\phi^{(0)}_{1}(\zeta)$ and $\phi^{(1)}_{1}(\zeta)$ can be extended to
a global entire function on the complex plane, when their derivatives and
themselves satisfy the condition

$G(E)\equiv\phi^{(0)}_{1}(\zeta;E)\phi^{(1)\prime}_{1}(\zeta;E)-\phi^{(1)}_{1}(\zeta;E)\phi^{(0)\prime}_{1}(\zeta;E),$
(45)

where $G(E)$ is the Wronskian of the local holomorphic series at the regular
singularities $0$ and $1$, regarded as a function of $E$, which vanishes on
quantized energies. This constitutes the primary result of this study.

###### Acknowledgements.

This study was supported by the National Natural Science Foundation of China
(Grant nos. 12104524).

##  Appendix A Bogoliubov transformation for the case $\Delta=0$

The energy spectrum for the quantum Rabi model with $\Delta=0$ can also be
obtained by using the Bogoliubov transformation. When $\Delta=0$, the
Hamiltonian is block diagonalized, where $H^{\prime}=\epsilon\sigma_{z}+\omega
a^{\dagger}a+g\sigma_{z}(a^{\dagger}+a)+\lambda\sigma_{z}(a^{\dagger
2}+a^{2})$. In such a case, $\sigma_{z}$ is a conserved quantity. For
$\sigma_{z}=1$, one obtains $H_{+}=\omega
a^{\dagger}a+g(a^{\dagger}+a)+\lambda(a^{\dagger 2}+a^{2})+\epsilon$, where
the subscript $\pm$ denotes the spin state with $\sigma_{z}=+1$ or $-1$
respectively. After performing a Bogoliubov transformation
$b_{+}=\mu_{+}a+\nu_{+}a^{\dagger}$ and
$b_{+}^{\dagger}=\nu_{+}^{*}a+\mu_{+}^{*}a^{\dagger}$ with
$|\mu_{+}|^{2}-|\nu_{+}|^{2}=1$, one obtains

$\displaystyle
H_{+}=\omega(|\mu|^{2}b^{\dagger}b+|\nu|^{2}bb^{\dagger})-\lambda(uv^{*}+u^{*}v)(b^{\dagger}b+bb^{\dagger})$
$\displaystyle+[(\lambda(\mu^{2}+\nu^{2})-\omega\mu\nu)b^{\dagger
2}+h.c.]+g[(\mu-\nu)b^{\dagger}+h.c.]+\epsilon.$ (46)

Hence, by choosing $\mu_{+}=\cosh r_{+}$, $\nu_{+}=\sinh r_{+}$, and $\tanh
2r_{+}=2\lambda/\omega$, one obtains

$H_{+}=\sqrt{\omega^{2}-4\lambda^{2}}\left(B_{+}^{\dagger}B_{+}+\frac{1}{2}\right)-\frac{g^{2}}{\omega+2\lambda}-\frac{\omega}{2}+\epsilon,$
(47)

where $B_{+}\equiv b_{+}+g(\omega-2\lambda)^{-1/4}(\omega+2\lambda)^{-3/4}$.
It implies that the positive branch spectrum of the generalized Rabi model
with $\Delta=0$ is given by

$E_{n_{+}}=\sqrt{\omega^{2}-4\lambda^{2}}\left(n_{+}+\frac{1}{2}\right)-\frac{g^{2}}{\omega+2\lambda}-\frac{\omega}{2}+\epsilon.$
(48)

Similarly, for $\sigma_{z}=-1$, one obtains $H_{-}=\omega
a^{\dagger}a-g(a^{\dagger}+a)-\lambda(a^{\dagger 2}+a^{2})-\epsilon$. Clearly,
by simply employing the transformation $g\rightarrow-g$,
$\lambda\rightarrow-\lambda$, and $\epsilon\rightarrow-\epsilon$, one would
obtain

$H_{-}=\sqrt{\omega^{2}-4\lambda^{2}}\left(B_{-}^{\dagger}B_{-}+\frac{1}{2}\right)-\frac{g^{2}}{\omega-2\lambda}-\frac{\omega}{2}-\epsilon,$
(49)

where $B_{-}\equiv b_{-}-g(\omega+2\lambda)^{-1/4}(\omega-2\lambda)^{-3/4}$.
Hence, the negative branch spectrum of the generalized Rabi model with
$\Delta=0$ is given by

$E_{n_{-}}=\sqrt{\omega^{2}-4\lambda^{2}}\left(n_{-}+\frac{1}{2}\right)-\frac{g^{2}}{\omega-2\lambda}-\frac{\omega}{2}-\epsilon.$
(50)

##  Appendix B Analytical approximation based on the second canonical form

After employing the transformations

$\displaystyle\phi_{1}(z)$
$\displaystyle=u_{1}(z)\exp\left\\{-\frac{1}{2}\bar{g}z-\frac{1}{4}\bar{\omega}z^{2}\right\\},$
(51a) $\displaystyle\phi_{2}(z)$
$\displaystyle=u_{2}(z)\exp\left\\{-\frac{1}{2}\bar{g}z+\frac{1}{4}\bar{\omega}z^{2}\right\\},$
(51b)

the coupled second-order differential equations becomes

$\displaystyle\frac{d^{2}u_{1}}{dz^{2}}+Q_{1}(z)u_{1}+\Delta_{1}(z)u_{2}$
$\displaystyle=0,$ (52a)
$\displaystyle\frac{d^{2}u_{2}}{dz^{2}}+Q_{2}(z)u_{2}+\Delta_{2}(z)u_{1}$
$\displaystyle=0.$ (52b)

where

$\displaystyle\Delta_{1}(z)\equiv\bar{\Delta}\exp\left(\frac{1}{2}\bar{\omega}z^{2}\right),\Delta_{2}(z)\equiv\bar{\Delta}\exp\left(-\frac{1}{2}\bar{\omega}z^{2}\right),$
(53a) $\displaystyle
Q_{1}(z)=\bar{\epsilon}-\frac{\bar{\omega}}{2}-\frac{\bar{g}^{2}}{4}-\bar{E}-\frac{1}{2}\bar{g}\left(\bar{\omega}-2\right)z-\frac{1}{4}\left(\bar{\omega}^{2}-4\right)z^{2},$
(53b) $\displaystyle
Q_{2}(z)=\bar{\epsilon}+\frac{\bar{\omega}}{2}-\frac{\bar{g}^{2}}{4}+\bar{E}+\frac{1}{2}\bar{g}\left(\bar{\omega}+2\right)z-\frac{1}{4}\left(\bar{\omega}^{2}-4\right)z^{2}.$
(53c)

The coupled second-order differential equations Eqs. (52a) and (52b) can be
reconfigured into a single fourth-order differential equation with respect to
$u_{1}$ or $u_{2}$

$\displaystyle\left(\frac{d^{2}}{dz^{2}}-2\bar{\omega}z\frac{d}{dz}+Q_{2}(z)+\bar{\omega}^{2}z^{2}-\bar{\omega}\right)\left(\frac{d^{2}}{dz^{2}}+Q_{1}(z)\right)u_{1}=\bar{\Delta}^{2}u_{1},$
(54a)
$\displaystyle\left(\frac{d^{2}}{dz^{2}}+2\bar{\omega}z\frac{d}{dz}+Q_{1}(z)+\bar{\omega}^{2}z^{2}+\bar{\omega}\right)\left(\frac{d^{2}}{dz^{2}}+Q_{2}(z)\right)u_{2}=\bar{\Delta}^{2}u_{2}.$
(54b)

Eq. (54a) can be reformulated as

$c_{0}(z)\frac{d^{4}u_{1}}{dz^{4}}+c_{1}(z)\frac{d^{3}u_{1}}{dz^{3}}+c_{2}(z)\frac{d^{2}u_{1}}{dz^{2}}+c_{3}(z)\frac{du_{1}}{dz}+c_{4}(z)u_{1}=0,$
(55)

where $c_{0}(z)\equiv 1$, $c_{1}(z)\equiv-2\bar{\omega}z$, $c_{2}(z)\equiv
Q_{1}+Q_{2}+\bar{\omega}^{2}z^{2}-\bar{\omega}$, $c_{3}(z)\equiv
2Q_{1}^{\prime}-2\bar{\omega}zQ_{1}$, and
$c_{4}(z)\equiv-\bar{\Delta}^{2}+Q_{1}^{\prime\prime}-2\bar{\omega}zQ_{1}^{\prime}+(Q_{2}+\bar{\omega}^{2}z^{2}-\bar{\omega})Q_{1}$.
A similar fourth-order differential equation can be derived from Eq. (54b)
through the transformation $\bar{\epsilon}\rightarrow-\bar{\epsilon}$,
$\bar{g}\rightarrow-\bar{g}$ and $\bar{\lambda}\rightarrow-\bar{\lambda}$. A
direct computation yields the following explicit expression of $c_{n}(z)$

$\displaystyle
c_{2}(z)=2\bar{\epsilon}-\bar{\omega}-\frac{\bar{g}^{2}}{2}+2\bar{g}z+\frac{1}{2}(\bar{\omega}^{2}+4)z^{2},$
(56a) $\displaystyle
c_{3}(z)=-\bar{g}(\bar{\omega}-2)+\left[-(\bar{\omega}^{2}-4)-2\bar{\omega}\bar{\epsilon}+\bar{\omega}^{2}+\frac{\bar{\omega}\bar{g}^{2}}{2}+2\bar{\omega}\bar{E}\right]z$
$\displaystyle+\bar{g}\bar{\omega}(\bar{\omega}-2)z^{2}+\frac{1}{2}\bar{\omega}(\bar{\omega}^{2}-4)z^{3},$
(56b) $\displaystyle
c_{4}(z)=-\frac{1}{2}(\bar{\omega}^{2}-4)+\left(\bar{\epsilon}-\frac{\bar{\omega}}{2}-\frac{\bar{g}^{2}}{4}\right)^{2}-\bar{E}^{2}-\bar{\Delta}^{2}$
$\displaystyle+\bar{g}\left(\bar{\omega}^{2}-(3+\bar{E})\bar{\omega}+2\bar{\epsilon}-\frac{\bar{g}^{2}}{2}\right)z$
$\displaystyle+\left[\left(\bar{\omega}-\frac{\bar{g}^{2}}{4}\right)\left(\bar{\omega}^{2}-4\right)+\frac{1}{2}\left(\bar{\epsilon}-\frac{\bar{\omega}}{2}-\frac{\bar{g}^{2}}{4}\right)\left(\frac{\bar{\omega}^{2}}{4}+1\right)-\bar{E}\bar{\omega}^{2}\right]z^{2}$
$\displaystyle-\frac{1}{2}\bar{g}(\bar{\omega}-2)(\bar{\omega}^{2}+\bar{\omega}+2)z^{3}-\frac{1}{16}(3\bar{\omega}^{2}+4)(\bar{\omega}^{2}-4)z^{4}.$
(56c)

Clearly, the polynomial $c_{n}(z)$ has the following properties: (1)
$c_{n}(z)$ represents an $n$-th order polynomials of $z$; (2)
$c_{n}(z)=O(\lambda^{-n})$. Hence, when $\lambda$ is a small but nonzero
parameter, one can keep the coefficients in the functions $c_{n}(z)$ up to
$\lambda^{-2}$, neglecting the two terms proportional to $c_{1}(z)$ and
$c_{0}(z)$, and obtaining the following second-order differential equation

$\displaystyle\frac{d^{2}u_{1}}{dz^{2}}+p_{1}(z)\frac{du_{1}}{dz}+q_{1}(z)u_{1}=0,$
(57a) $\displaystyle
p_{1}(z)\equiv\frac{c_{3}(z)}{c_{2}(z)},q_{1}(z)\equiv\frac{c_{4}(z)}{c_{2}(z)},$
(57b)

where

$\displaystyle c_{2}(z)$
$\displaystyle\approx-\frac{\bar{g}^{2}}{2}+\frac{\bar{\omega}^{2}}{2}z^{2}=\frac{\bar{\omega}^{2}}{2}\left(z+\frac{\bar{g}}{\bar{\omega}}\right)\left(z-\frac{\bar{g}}{\bar{\omega}}\right),$
$\displaystyle c_{3}(z)$
$\displaystyle\approx-\bar{\omega}\bar{g}+\bar{\omega}\left(2\bar{E}-2\bar{\epsilon}+\frac{\bar{g}^{2}}{2}\right)z+\bar{g}\bar{\omega}(\bar{\omega}-2)z^{2}+\frac{1}{2}\bar{\omega}^{3}z^{3},$
$\displaystyle c_{4}(z)$
$\displaystyle\approx-\frac{1}{2}\bar{\omega}^{2}+\left(\bar{\epsilon}-\frac{\bar{\omega}}{2}-\frac{\bar{g}^{2}}{4}\right)^{2}-\bar{E}^{2}-\bar{\Delta}^{2}$
$\displaystyle+\bar{g}\left(\bar{\omega}^{2}-(3+\bar{E})\bar{\omega}+2\bar{\epsilon}-\frac{\bar{g}^{2}}{2}\right)z$
$\displaystyle+\left\\{\bar{\omega}^{3}+\left[\frac{1}{8}\left(\bar{\epsilon}-\frac{\bar{\omega}}{2}-\frac{\bar{g}^{2}}{4}\right)-\left(\frac{g^{2}}{4}+\bar{E}\right)\right]\bar{\omega}^{2}+\frac{7}{8}\bar{g}^{2}\right\\}z^{2}$
$\displaystyle-\frac{1}{2}\bar{g}\bar{\omega}^{2}(\bar{\omega}-1)z^{3}-\frac{1}{16}\bar{\omega}^{2}(3\bar{\omega}^{2}-8)z^{4}.$

A direct computation yields the partial fraction decomposition of $p_{1}(z)$
and $q_{1}(z)$

$\displaystyle p_{1}(z)$
$\displaystyle=\alpha_{1}z+\alpha_{2}+\frac{\beta_{1}}{z-\bar{g}/\bar{\omega}}+\frac{\beta_{2}}{z+\bar{g}/\bar{\omega}},$
(58a) $\displaystyle q_{1}(z)$
$\displaystyle=\gamma_{1}z^{2}+\gamma_{2}z+\gamma_{3}+\frac{\delta_{1}}{z-\bar{g}/\bar{\omega}}+\frac{\delta_{2}}{z+\bar{g}/\bar{\omega}}.$
(58b)

where

$\displaystyle\alpha_{1}$ $\displaystyle\equiv\bar{\omega},\alpha_{2}\equiv
2\bar{g}\left(1-\frac{2}{\bar{\omega}}\right),\gamma_{1}\equiv
1-\frac{3\bar{\omega}^{2}}{8},\gamma_{2}\equiv-\bar{g}(\bar{\omega}-1),$ (59a)
$\displaystyle\beta_{1}$
$\displaystyle\equiv\frac{2}{\bar{\omega}}\left[\bar{E}-\bar{\epsilon}+\bar{g}^{2}\left(1-\frac{1}{\bar{\omega}}\right)\right]-1,\beta_{2}\equiv\frac{2}{\bar{\omega}}\left[\bar{E}-\bar{\epsilon}+\frac{\bar{g}^{2}}{\bar{\omega}}\right]+1,$
(59b) $\displaystyle\gamma_{3}$
$\displaystyle\equiv\frac{11}{4}\frac{\bar{g}^{2}}{\bar{\omega}^{2}}+\frac{15}{8}\bar{\omega}-2\bar{E}+\frac{\bar{\epsilon}}{4}-\frac{15}{16}\bar{g}^{2},$
(59c) $\displaystyle\delta_{1}$
$\displaystyle\equiv\frac{\bar{g}}{\bar{\omega}}\left(\frac{31}{16}\bar{\omega}+\frac{\bar{\epsilon}}{8}-3-2\bar{E}+\frac{2\bar{\epsilon}}{\bar{\omega}}\right)+\frac{\bar{g}^{3}}{\bar{\omega}^{3}}\left(\frac{13}{32}\bar{\omega}^{2}+\frac{11}{8}\right)$
$\displaystyle-\frac{\bar{\omega}^{2}}{2}+\left(\bar{\epsilon}-\frac{\bar{\omega}}{2}-\frac{\bar{g}^{2}}{4}\right)^{2}-\bar{E}^{2}-\bar{\Delta}^{2},$
(59d) $\displaystyle\delta_{2}$
$\displaystyle\equiv\frac{\bar{g}}{\bar{\omega}}\left(\frac{1}{16}\bar{\omega}-\frac{\bar{\epsilon}}{8}+3+2\bar{E}-\frac{2\bar{\epsilon}}{\bar{\omega}}\right)-\frac{\bar{g}^{3}}{\bar{\omega}^{3}}\left(\frac{1}{32}\bar{\omega}^{2}+\frac{11}{8}\right)$
$\displaystyle+\frac{\bar{\omega}^{2}}{2}-\left(\bar{\epsilon}-\frac{\bar{\omega}}{2}-\frac{\bar{g}^{2}}{4}\right)^{2}+\bar{E}^{2}+\bar{\Delta}^{2}.$
(59e)

After the transformation

$u_{1}(z)\equiv
U_{1}(z)\left(z-\frac{\bar{g}}{\bar{\omega}}\right)^{-\beta_{1}/2}\left(z+\frac{\bar{g}}{\bar{\omega}}\right)^{-\beta_{2}/2}\exp\left\\{-\frac{1}{4}\alpha_{1}z^{2}-\frac{1}{2}\alpha_{2}z\right\\},$
(60)

One derives the secondary form of the differential equation Eq. (57a), namely

$\displaystyle
U_{1}^{\prime\prime}=-\left(\lambda_{1}z^{2}+\lambda_{2}z+\lambda_{3}+\frac{\mu_{1}}{z-q}+\frac{\mu_{2}}{z+q}\right.$
$\displaystyle\left.+\frac{\nu_{1}}{(z-q)^{2}}+\frac{\nu_{2}}{(z+q)^{2}}\right)U_{1},$
(61)

where $q\equiv\bar{g}/\bar{\omega}$ and

$\displaystyle\lambda_{1}\equiv\gamma_{1}-\frac{\alpha_{1}}{4},\lambda_{2}\equiv\gamma_{2}-\frac{1}{2}\alpha_{1}\alpha_{2},$
$\displaystyle\lambda_{3}\equiv\gamma_{3}-\frac{\alpha_{1}}{2}\left(1+\beta_{1}+\beta_{2}\right)-\frac{\alpha_{2}^{2}}{4},$
$\displaystyle\mu_{1}\equiv\delta_{1}-\frac{1}{2}\left(\frac{\bar{g}}{\bar{\omega}}\alpha_{1}+\alpha_{2}+\frac{\bar{\omega}}{4\bar{g}}\beta_{2}\right)\beta_{1},$
$\displaystyle\mu_{2}\equiv\delta_{2}+\frac{1}{2}\left(\frac{\bar{g}}{\bar{\omega}}\alpha_{1}-\alpha_{2}+\frac{\bar{\omega}}{4\bar{g}}\beta_{1}\right)\beta_{2},$
$\displaystyle\nu_{1}\equiv\frac{\beta_{1}}{2}\left(1-\frac{\beta_{1}}{2}\right),\nu_{2}\equiv\frac{\beta_{2}}{2}\left(1-\frac{\beta_{2}}{2}\right).$
(62)

Notably, for the special case where $g=0$, one obtains

$\frac{d^{2}U_{1}}{dz^{2}}+\left(\lambda_{1}z^{2}+\lambda_{2}z+\lambda_{3}+\frac{\mu}{z}+\frac{\nu}{z^{2}}\right)U_{1}=0,$
(63)

where

$\displaystyle\lambda_{1}=\left(1+\frac{\bar{\omega}}{2}\right)\left(1-\frac{3\bar{\omega}}{4}\right),\lambda_{2}=0,\lambda_{3}=\frac{11}{8}\bar{\omega}-4\bar{E}+\frac{9\bar{\epsilon}}{4},$
$\displaystyle\mu=0,\nu\equiv\frac{2}{\bar{\omega}}(\bar{E}-\bar{\epsilon})\left[1-\frac{2}{\bar{\omega}}(\bar{E}-\bar{\epsilon})\right].$
(64)

Intriguingly, Eq. (63) can be expressed in the the second normal form of the
biconfluent Heun equation as

$U_{1}^{\prime\prime}=\left(\frac{\zeta^{2}}{4}+\frac{\delta}{2}\zeta-\frac{1+2\alpha-\gamma}{2}+\frac{\delta^{2}}{4}+\frac{\gamma\delta+\beta}{2\zeta}+\frac{\gamma(2+\gamma)}{4\zeta^{2}}\right)U_{1},$
(65)

where $\xi\equiv e^{-i\pi/4}(4\lambda_{1})^{1/4}z$, $\beta=\delta=0$, and

$\alpha\equiv-2\left(\frac{1}{\bar{\omega}}+2\right)\bar{E}+\left(\frac{9}{4}+\frac{2}{\bar{\omega}}\right)\bar{\epsilon}+\frac{11}{8}\bar{\omega}-\frac{1}{2},\gamma\equiv-\frac{4}{\bar{\omega}}(\bar{E}-\bar{\epsilon}).$
(66)

After the transformation

$V_{1}(\zeta)\equiv
U_{1}(\zeta)\zeta^{\gamma/2}\exp\left(\frac{\delta}{2}\zeta+\frac{1}{4}\zeta^{2}\right),$
(67)

one obtains the first normal form of the biconfluent Heun equation

$\frac{d^{2}V_{1}}{d\zeta^{2}}-\left(\frac{\gamma}{\zeta}+\delta+\zeta\right)\frac{dV_{1}}{d\zeta}+\left(\alpha-\frac{\beta}{\zeta}\right)V_{1}=0.$
(68)

Eq. (68) can be solved by the biconfluent Heun function
$\mbox{BCH}(\alpha,\beta,\gamma,\delta;\zeta)$, which has the following power
series

$\mbox{BCH}(\alpha,\beta,\gamma,\delta;\zeta)=\sum_{n=0}^{\infty}\frac{A_{n}}{(-\gamma)^{(n)}}\frac{\zeta^{n}}{n!},$
(69)

where $x^{(n)}\equiv x(x+1)\cdots(x+n-1)$ is the rising factorial which
satisfies $x^{(n+1)}=(x+n)x^{(n)}$ and $x^{(0)}\equiv 1$, and the coefficients
$A_{n}$ for $n\geq 2$ are determined by the three-term recurrence relation

$A_{n+2}-(\delta(n+1)+\beta)A_{n+1}+(n+1)(n-\gamma)(\alpha-n)A_{n},$ (70)

with $A_{0}=1$ and $A_{1}=\beta$. For the special case where $g=0$, one
obtains $\beta=0$ and $\delta=0$, and thus the recurrence relation is
simplified to

$A_{n+2}=(n+1)(\gamma-n)(\alpha-n)A_{n},$ (71)

## References

  * (1) I. I. Rabi, Phys. Rev. 49, 324 (1936).
  * (2) I. I. Rabi, Phys. Rev. 51, 652 (1937).
  * (3) S. Gleyzes, S. Kuhr, C. Guerlin, J. Bernu, S. Deleglise, U. Busk Hoff, M. Brune, J. M. Raimond, and S. Haroche, Nature 446, 297 (2007).
  * (4) C. Guerlin, B. Julien, S. Deleglise, C. Sayrin, S. Gleyzes, S. Kuhr, M. Brune, J. M. Raimond, and S. Haroche, Nature 448, 889 (2007).
  * (5) P. Forn-Díaz, J. Lisenfeld, D. Marcos, J. J. Garcia-Ripoll, E. Solano, C. J. P. M. Harmans, and J. E. Mooij, Phys. Rev. Lett. 105, 237001 (2010).
  * (6) M. D. LaHaye, J. Suh, P. M. Echternach, K. C. Schwab, and M. L. Roukes, Nature 459, 960 (2009).
  * (7) L. Yanbing, and A. A. Houck, Nat. Phys. 13, 48 (2017).
  * (8) D. Braak, Phys. Rev. Lett. 107, 100401 (2011).
  * (9) V. Bargmann, Commun. Pure Appl. Math. 14, 187 (1961).
  * (10) I. E. Segal, Mathematical problems of relativistic physics (American Mathematical Soc., 1963).
  * (11) W. M. Zhang, R. Gilmore, and D. H. Feng, Rev. Mod. Phys. 4, 867 (1990).
  * (12) C. F. Kam, W. M. Zhang, and D. H. Feng, Coherent States: New Insights into Quantum Mechanics with Applications (Springer Nature, 2023).
  * (13) H. H. Zhong, Q. T. Xie, M. T. Batchelor, and C. H. Lee, J. Phys. A Math. Theor. 46, 415302 (2013).
  * (14) A. Ronveaux and F. M. Arscott, Heun’s differential equations (Oxford University Press, 1995).
  * (15) B. R. Judd, J. Phys. C: Solid State Phys. 12, 1685 (1979).
  * (16) W. Pauli, Pauli lectures on physics. V. Wave mechanics (Dover Publications, 2000).
  * (17) A. M. Ishkhanyan, and A. E. Grigoryan, J. Phys. A Math. Theor. 47, 465205 (2014).
  * (18) A. M. Ishkhanyan, T. A. Shahverdyan, and T. A. Ishkhanyan, Eur. Phys. J. D. 69, 1 (2015).
  * (19) C. F. Kam, and Y. Chen, New J. Phys. 22, 023021 (2020).
  * (20) C. F. Kam, and Y. Chen, Z. Angew. Math. Phys. 72, 1 (2021).
  * (21) S. Dorosz, T. Platini, and D. Karevski, Phys. Rev. E 77, 051120 (2008).
  * (22) C. A. Downing, and M. E. Portnoi, Phys. Rev. B 94, 165407 (2016).
  * (23) C. F. Kam, and Y. Chen, Ann. Phys. (Berl.) 533, 2000349 (2021).
  * (24) C. F. Kam, and Y. Chen, arXiv preprint arXiv:2301.04816
  * (25) C. C. Gerry, Phys. Rev. A 37, 2683 (1988).
  * (26) V. Del Valle, S. Zippilli, F. P. Laussy, A. Gonzalez-Tudela, G. Morigi, and C. Tejedor, Phys. Rev. B 81, 035302 (2010).
  * (27) M. Tse, H. Yu, N. Kijbunchoo, A. Fernandez-Galiana, P. Dupej, L. Barsotti, C. D. Blair, D. D. Brown, S. E. Dwyer, A. Effler, and others, Phys. Rev. Lett. 123, 231107 (2019).
  * (28) F. Acernese, M. Agathos, L. Aiello, A. Allocca, A. Amato, S. Ansoldi, S. Antier, M. Arène, N. Arnaud, S. Ascenzi, and others, Phys. Rev. Lett. 123, 231108 (2019).
  * (29) S. E. Dwyer, G. L. Mansell, and L. McCuller, Galaxies 10, 46 (2022).
  * (30) T. Gehring, V. Händchen, J. Duhme, F. Furrer, T. Franz, C. Pacher, R. F. Werner, and R. Schnabel, Nat. Commun. 6, 1 (2015).
  * (31) B. J. Lawrie, P. D. Lett, A. M. Marino, and R. C. Pooser, ACS photonics 6, 1307 (2019).
  * (32) J. Junker, D. Wilken, E. Huntington, and M. Heurs, Opt. Express 29, 6053 (2021).
  * (33) C. F. Kam, and X. D. Hu, arXiv preprint arXiv:2312.10820
  * (34) C. F. Kam, and X. D. Hu, arXiv preprint arXiv:2401.03617