Analytic solutions of 2-photon and two-mode Rabi models
{Yao-Zhong Zhang}$\
we obtain analytic representations of solutions to the 2-photon and two-mode quantum Rabi models. In each case
a transcendental function is analytically derived whose zeros give the energy spectrum of the model. The zeros can be numerically found by standard root-search techniques. We also present analytic solution to the driven Rabi model with broken ${\cal Z}_2$ symmetry. \end{abstract} \vskip.1in {PACS numbers}: 03.65.Ge
02.30.Ik
42.50.Pq % 03.65Fd - Algebraic methods %{Keywords}: Quasi-exactly solvable systems
Rabi model and generalizations
%\vspace{0.5cm} %\end{titlepage} %02.30.lk Integrable systems %03.65.Fd - Algebraic methods %03.65.Ge Solutions of wave equations: bound state %03.75.-b Matter waves (for atoms) %03.75.Kk Dynamical properties of condensates collective and hydrodynamic excitations
superfluid flow %03.75.Mn Multicompoment condensates
spinor condensates %05.30.Jp Boson systems %37.10.Jk Atoms in optical lattice %37.10.Gh Atoms traps and guides %42.50.-p - Quantum optics %42.50.Pq - Cavity quantum electrodynamics; micromasers %67.85.-d Ultracold gases
trapped gases %67.85.De Dynamical properties of condensates
excitations
superfluid flow %67.85.Fg Multicomponent condensates
spinor condesates %67.85.Hj Bose-Einstein condensation in optical lattice %\setcounter{section}{0} %\setcounter{equation}{0} \section{Introduction} Spin-boson systems describe interactions between spins and harmonic oscillators (boson modes) and have played a prominent role in modeling the ubiquitous matter-light interactions in modern physics. One of the well-known spin-boson systems is the quantum Rabi model
introduced by Rabi over 70 years ago. This model describes the interaction of a two-level atom with a harmonic mode of quantized electromagnetic field. Due to the simplicity of its Hamiltonian
the Rabi model has served as the basis for understanding matter-light interactions
and has a variety of applications ranging from quantum optics \cite{Vedral06} to solid state semiconductor systems \cite{Khitrova06} and molecular physics \cite{Thanopulos04}. It has also played a significant role in the novel research field of cavity and circuit quantum electrodynamics \cite{Englund07,Niemczyk10}. The Hamiltonian of the Rabi model is given by \beq H_R=\omega b^\dagger b+\Delta\
\sigma_z+g\
\sigma_x(b^\dagger+b)
\eeq where $g$ is the spin-boson interaction strength
$\sigma_z
\sigma_x$ are the Pauli matrices describing the two atomic levels separated by energy difference $2\Delta$
and $b^\dagger$ ($b$) are creation (annihilation) operators of a boson mode with frequency $\omega$. In the Bargmann space with the monomials $\left\{\frac{w^n}{\sqrt{n!}}\right\}$ as basis vectors
the boson creation and annihilation operators can be realized as $b^\dagger = w$
$b=\frac{d}{dw}$. In terms of the two-component wavefunction $\psi(w)=\left(\psi_+(w)
\psi_-(w)\right)^T$
the time-independent Schr\"odinger equation of the model % $H_R\psi(w) =E\psi(w)$
gives a system of two coupled linear differential equations with rational coefficients. The system depends on a spectral parameter $E$. The wavefunction components are elements of Bargmann-Hilbert space of entire functions of one variable $w\in {C}$ \cite{Schweber67}. %harmonic creation and annihilation operators $a^\dagger
a$ %are realized as $a^\dagger\rightarrow z
~a\rightarrow \frac{d}{dz}$. The scalar product of any two elements $f(w)
g(w)$ in the Bargmann-Hilbert space is given by \beq (f
g)=\int\
\overline{f(w)}\
g(w)\
d\mu(w)
\label{inner-product1} \eeq where $d\mu(w)=\frac{1}{\pi}e^{-|w|^2}\
dx\
dy$. %An entire function belongs to the Bargmann-Hilbert space if it has proper growth at infinity. The energy $E$ belongs to the spectrum of the problem if and only if for this value of $E$ the system of coupled differential equations (given by the time-independent Schr\"odinger equation) has entire solution \cite{Schweber67}. By means of the Bargmann-Hilbert space approach
Braak \cite{Braak11} recently presented a pair of transcendental functions defined as infinite power series expansions with coefficients satisfying three-term recurrence relations
and argued that the spectrum of the Rabi model is given by the zeros of the transcendental functions. This theoretical progress has renewed the interest in the Rabi and related models \cite{Solano11}-\cite{Moroz14}. In this paper we solve the 2-photon and two-mode generalizations of the quantum Rabi model. The 2-photon and two-mode Rabi models are phenomenological models describing a two-level atom interacting with 2 photons and 2 harmonic modes
respectively. They can be experimentally realized in circuit quantum electromagnetic systems \cite{Niemczyk10}
and have established applications in many research fields
including Rubidium atoms \cite{Brune87} and quantum dots \cite{Valle10,Ota11}. The reset of this paper is as follows. Applying the Bogoliubov type transformations and Bargmann-Hilbert spaces
we obtain analytic representations of solutions to the two models in sections \ref{2-photon} and \ref{two-mode}
respectively. For each case
a transcendental function is analytically found whose zeros give the energy eigenvalues of the model for arbitrary coupling strength and detuning. We draw conclusions in section \ref{summary}. %The zeros can be found numerically by standard root-search techniques %(see e.g. \cite{Leaver86,Liu92} and references therein). In the Appendix
we apply our approach to the driven Rabi model without ${\cal Z}_2$ symmetry and derive a continued fraction equation for its energy spectrum. %\setcounter{section}{1} %\setcounter{equation}{0} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{2-photon quantum Rabi model}\label{2-photon} The Hamiltonian of the 2-photon Rabi model reads \begin{equation} H_{2p}=\omega b^\dagger b+\Delta\sigma_z+g\
\sigma_x\left[(b^\dagger)^2 +b^2\right].\label{2-photon-RabiH} \end{equation} %where $g$ is the interaction strength
%$\sigma_z
\sigma_x$ are the Pauli matrices describing the two atomic levels separated %by energy difference $2\Delta$
and $b^\dagger$ ($b$) are %creation (annihilation) operators of a boson mode with frequency $\omega$. Let us make a canonical Bogoliubov transformation from $b
b^\dagger$ to the squeezed bosons $a
a^\dagger$ \cite{Emary02}
\beq b=\frac{a+\tau\, a^\dagger}{\sqrt{1-\tau^2}}
~~~~~~b^\dagger=\frac{\tau\,a+a^\dagger}{\sqrt{1-\tau^2}}
\eeq where $|\tau|<1$ is a real parameter. In terms of the squeezed bosons
the Hamiltonian (\ref{2-photon-RabiH}) takes the form \beqa \tilde{H}_{2p}&=&\Delta\sigma_z+\frac{1}{1-\tau^2}\left[\left(\omega\tau+g\sigma_x(1+\tau^2)\right)\left((a^\dagger)^2 +a^2\right)\right.\n & &\left.+\left(\omega(1+\tau^2)+4g\tau\sigma_x\right)a^\dagger a+\omega\tau^2+2g\tau\sigma_x\right]. \label{2-photon-RabiH1} \eeqa Introduce the operators $K_\pm
K_0$ \beq K_+=\frac{1}{2}(a^\dagger)^2
~~~~K_-=\frac{1}{2}a^2
~~~~K_0=\frac{1}{2}\left(a^\dagger a+\frac{1}{2}\right). \label{2-photon-Rabi-boson} \eeq Then (\ref{2-photon-RabiH1}) becomes %\begin{widetext} %\beq %\tilde{H}_{2-p}=\left( %\begin{array}{cc} %2\frac{\omega\tau+g(1+\tau^2)}{1-\tau^2} (K_++K_-) & \Delta
%+2\frac{\omega(1+\tau^2)+4g\tau}{1-\tau^2}K_0-\frac{1}{2}\omega
Applying Bogoliubov transformations and Bargmann-Hilbert spaces, we obtain analytic representations of solutions to the 2-photon and two-mode quantum Rabi models. In each case, a transcendental function is analytically derived whose zeros give the energy spectrum of the model. The zeros can be numerically found by standard root-search techniques. We also present analytic solution to the driven Rabi model with broken ${\cal Z}_2$ symmetry.
1304.7827.tex

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Latex twice to get the right references.                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentstyle[12pt]{article}

\renewcommand{\baselinestretch}{1.2} % to make document double-spaced
%\oddsidemargin0mm \evensidemargin-4.5mm \topmargin-2.0cm %-27mm
\oddsidemargin14mm \evensidemargin14mm \topmargin0mm % Gustav's offsets
%\pagestyle{myheadings}
\textwidth=160mm
\textheight=240mm % 250mm
\hoffset -15mm
\voffset -20mm

\renewcommand{\thefootnote}{\fnsymbol{footnote}}

%\textwidth=175mm
%\textheight=260mm
%\hoffset -20mm
%\voffset -55mm


%\textwidth=165mm
%\textheight=250mm
%\hoffset -15mm
%\voffset -50mm

% \pagestyle{myheadings}
%\markright{~~~~~~~~~~~ Zhang et al: Primary fields of $gl(2|2)$
%conformal field theory}


\makeatletter



%%%%%%%%%%%%%%Equation version%%%%%%%%%%

\newcommand {\beq}{\begin{equation}}
\newcommand {\eeq}{\end{equation}}
\newcommand {\beqa}{\begin{eqnarray}}
\newcommand {\eeqa}{\end{eqnarray}}         %Equation version
\newcommand {\beqs}{\begin{eqnarray*}}
\newcommand {\eeqs}{\end{eqnarray*}}
\newcommand {\bds}{\begin{displaymath}}
\newcommand {\eds}{\end{displaymath}}
\newcommand {\n}{\nonumber\\}
\newcommand {\nn}{\nonumber}
%%%%%%%%%%%%%%%

%%%%%%%%%%%%%Reference version%%%%%%%%%%%%%%%%
\newcommand {\bebb}{\begin{thebibliography}}
\newcommand {\eebb}{\end{thebibliography}}      %Reference version
\newcommand {\bbit}{\bibitem}
%\newcommand {\Bibitem}[1]{\bibitem{#1}}
%\newcommand {\Bibitem}[1]{\bibitem{#1} {\tt #1}}

\def\dz{\frac{d}{dz}}

\newcommand{\cA}{{\cal A}}
\newcommand{\cB}{{\cal B}}
\newcommand{\cC}{{\cal C}}
\newcommand{\cD}{{\cal D}}
\newcommand{\cF}{{\cal F}}
\newcommand{\cG}{{\cal G}}
\newcommand{\cH}{{\cal H}}
\newcommand{\cI}{{\cal I}}
\newcommand{\cK}{{\cal K}}
\newcommand{\cM}{{\cal M}}
\newcommand{\cN}{{\cal N}}
\newcommand{\cP}{{\cal P}}
\newcommand{\cQ}{{\cal Q}}
\newcommand{\cR}{{\cal R}}
\newcommand{\cS}{{\cal S}}
\newcommand{\cT}{{\cal T}}
\newcommand{\cU}{{\cal U}}
\newcommand{\cV}{{\cal V}}
\newcommand{\cW}{{\cal W}}
\newcommand{\cY}{{\cal Y}}
\newcommand{\cZ}{{\cal Z}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\evn}{{\overline{0}}}
\newcommand{\odd}{{\overline{1}}}

%\newcommand
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Useful symbols         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\a{\alpha}
\def\bt{\beta}

\def\d{\delta}
\def\D{\Delta}

\def\eps{\epsilon}
\def\veps{\varepsilon}


\def\G{\Gamma}
\def\gm{\gamma}

\def\k{\kappa}

\def\l{\lambda}
\def\L{\Lambda}

\def\o{\omega}
\def\O{\Omega}

\def\vph{\varphi}

\def\vpi{\varpi}

\def\s{\sigma}
\def\S{\Sigma}
\def\vsgm{\varsigma}

\def\t{\theta}
\def\T{\Theta}
\def\vt{\vartheta}

\def\z{\zeta}


%%%%%%%%%%%%%%%For over letters%%%%%%%%%%%%%%%%%%%%
\def\tl{\tilde}  %%%%%%%%%%For over letters%%%%%%%%%
\def\wht{\widehat}
\def\arr{\arrow}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\empst{\emptyset}
\def\prm{\prime}
\def\exst{\exists}
\def\el{\ell}
\def\nba{\nabla}
\def\p{\partial}
\def\bksh{\backslash}
\def\sd{\surd}
\def\nft{\infty}
\def\fal{\forall}
\def\cpd{\coprod}
\def\pd{\prod}
\def\bgp{\bigup}
\def\bgc{\bigcap}
\def\prt{\propto}

%%%%%%%%%%%%%%%%%%%%brackets
\def\la{\langle}
\def\ra{\rangle}
\def\rb{\rbrace}
\def\lb{\lbrace}
\def\lt{\left}
\def\rt{\right}
%%%%%%%%%%%%%%%%%%%%%%sinh
\def\sh{\sinh}
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%\ldots low %%%\vdots vertical
%%%%%%\cdots center%%%%%ddots diagonal


\def\tri{\triangle}
\def\ot{\otimes}
\def\bot{\bigotimes}
\def\op{\oplus}
\def\bop{\bigoplus}

\def\rtarr{\rightarrow}
\def\Rtarr{\Rightarrow}
\def\Ltarr{\Leftarrow}
\def\ltarr{\leftarrow}
\def\lrtarr{\leftrightarrow}
\def\mpto{\mapsto}

\def\dg{\dagger}
\def\apx{\approx}
\def\eqv{\equiv}
\def\sbs{\sebset}
\def\sbsq{\sebseteq}
\def\smq{\simeq}

\def\psd{\psi ^{\dagger}}
\def\Psd{\Psi ^{\dagger}}
\def\psb{\bar{\psi}}
\def\psbd{\bar{\psi} ^{\dagger}}
\def\psp{\psi ^{\prime}}
\def\pspd{{\psi ^{\prime}}^{\dagger}}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    End of Yao's personal definition                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% My definitions and modification
\def\dz{\frac{d}{dz}}

\begin{document}


%\begin{titlepage}

\begin{flushright}
\end{flushright}

%\baselineskip =16pt

\vskip 1cm

\begin{center}
%\title
{\large\bf Analytic solutions of 2-photon and two-mode Rabi models} % in Bargmann-Hilbert spaces

\vspace{1cm}

%\author{
{Yao-Zhong Zhang}$\,$\footnote{Email: yzz@maths.uq.edu.au}
\vskip.1in

{\em School of Mathematics and Physics,
The University of Queensland, Brisbane, Qld 4072, Australia}


\end{center}

\date{}

%\maketitle

%\vspace{2cm}

\begin{abstract}
Applying Bogoliubov transformations and Bargmann-Hilbert spaces, we obtain analytic representations of 
solutions to the 2-photon and two-mode quantum Rabi models. In each case, a transcendental function is
analytically derived whose zeros give the energy spectrum of the model. The zeros can be numerically 
found by standard root-search techniques. We also present analytic solution to the driven Rabi model
with broken ${\cal Z}_2$ symmetry.  
\end{abstract}

\vskip.1in

{\it PACS numbers}:  03.65.Ge, 02.30.Ik, 42.50.Pq
% 03.65Fd - Algebraic methods

%{\it Keywords}: Quasi-exactly solvable systems, Rabi model and generalizations,

%\vspace{0.5cm}

%\end{titlepage}


%02.30.lk Integrable systems

%03.65.Fd - Algebraic methods
%03.65.Ge Solutions of wave equations: bound state
%03.75.-b Matter waves (for atoms)
%03.75.Kk Dynamical properties of condensates collective and hydrodynamic excitations, superfluid flow
%03.75.Mn Multicompoment condensates, spinor condensates

%05.30.Jp Boson systems

%37.10.Jk Atoms in optical lattice
%37.10.Gh Atoms traps and guides

%42.50.-p - Quantum optics
%42.50.Pq - Cavity quantum electrodynamics; micromasers 

%67.85.-d Ultracold gases, trapped gases
%67.85.De Dynamical properties of condensates, excitations, superfluid flow
%67.85.Fg Multicomponent condensates, spinor condesates
%67.85.Hj Bose-Einstein condensation in optical lattice




%\setcounter{section}{0}
%\setcounter{equation}{0}

\section{Introduction}
Spin-boson systems describe interactions between spins and harmonic oscillators (boson modes)
and have played a prominent role in modeling the ubiquitous matter-light interactions in modern physics.
One of the well-known spin-boson systems is the quantum Rabi model,
introduced by Rabi over 70 years ago.  This model describes the interaction
of a two-level atom with a harmonic mode of quantized electromagnetic field. 
Due to the simplicity of its Hamiltonian, the Rabi model has served as the basis
for understanding matter-light interactions, and has a variety of applications ranging from  
quantum optics \cite{Vedral06} to solid state semiconductor systems \cite{Khitrova06} 
and molecular physics \cite{Thanopulos04}. It has also played a significant role in
the novel research field of cavity and circuit quantum electrodynamics \cite{Englund07,Niemczyk10}. 

The Hamiltonian of the Rabi model is given by
\beq
H_R=\omega b^\dagger b+\Delta\,\sigma_z+g\,\sigma_x(b^\dagger+b),
\eeq
where $g$ is the spin-boson interaction strength, 
$\sigma_z, \sigma_x$ are the Pauli matrices describing the two atomic levels separated 
by energy difference $2\Delta$, and $b^\dagger$ ($b$) are 
creation (annihilation) operators of a boson mode with frequency $\omega$. 
In the Bargmann space with the monomials $\left\{\frac{w^n}{\sqrt{n!}}\right\}$ as basis 
vectors, the boson creation and annihilation operators can be realized as 
$b^\dagger = w$, $b=\frac{d}{dw}$.
In terms of the two-component wavefunction $\psi(w)=\left(\psi_+(w), \psi_-(w)\right)^T$,
the time-independent Schr\"odinger equation of the model  % $H_R\psi(w) =E\psi(w)$,
gives a system of two coupled linear differential equations with rational coefficients. 
The system depends on a spectral parameter $E$. The wavefunction components are
elements of Bargmann-Hilbert space of entire functions of one variable $w\in {\bf C}$ \cite{Schweber67}.
%harmonic creation and annihilation operators $a^\dagger, a$ 
%are realized as $a^\dagger\rightarrow z,~a\rightarrow \frac{d}{dz}$.
The scalar product of any two elements $f(w), g(w)$ in the Bargmann-Hilbert space is given by
\beq
(f,g)=\int\,\overline{f(w)}\,g(w)\,d\mu(w),\label{inner-product1}
\eeq
where $d\mu(w)=\frac{1}{\pi}e^{-|w|^2}\,dx\,dy$. 
%An entire function belongs to the Bargmann-Hilbert space if it has proper growth at infinity. 
The energy $E$ belongs to the spectrum
of the problem if and only if for this value of $E$ the system of coupled differential equations
(given by the time-independent Schr\"odinger equation) has entire solution \cite{Schweber67}.

By means of the Bargmann-Hilbert space approach, Braak \cite{Braak11} recently presented a pair of 
transcendental functions defined as infinite 
power series expansions with coefficients satisfying three-term recurrence relations, 
and argued that the spectrum of the Rabi model is given by the zeros of the transcendental functions. 
This theoretical progress has renewed the interest in the Rabi and related models
\cite{Solano11}-\cite{Moroz14}.

In this paper we solve the 2-photon and two-mode generalizations of the quantum Rabi model. 
The 2-photon and two-mode Rabi models
are phenomenological models describing a two-level atom interacting with 2 photons and 2 harmonic modes,
respectively. They can be experimentally realized in circuit quantum electromagnetic systems
\cite{Niemczyk10}, and have established applications in many research fields,
including Rubidium atoms \cite{Brune87} and  quantum dots \cite{Valle10,Ota11}. 

The reset of this paper is as follows. 
Applying the Bogoliubov type transformations and Bargmann-Hilbert spaces, we 
obtain analytic representations of solutions to the two models in sections \ref{2-photon} and
\ref{two-mode}, respectively. For each case, a transcendental function is
analytically found whose zeros give the energy eigenvalues of the model for arbitrary 
coupling strength and detuning. We draw conclusions in section \ref{summary}. 
%The zeros can be found numerically by standard root-search techniques
%(see e.g. \cite{Leaver86,Liu92} and references therein). 
In the Appendix, we apply our approach to the driven Rabi model without ${\cal Z}_2$ symmetry 
and derive a continued fraction equation for its energy spectrum.

  



%\setcounter{section}{1}
%\setcounter{equation}{0}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{2-photon quantum Rabi model}\label{2-photon}
The Hamiltonian of the 2-photon Rabi model reads
\begin{equation}
H_{2p}=\omega b^\dagger b+\Delta\sigma_z+g\,\sigma_x\left[(b^\dagger)^2 +b^2\right].\label{2-photon-RabiH}
\end{equation}
%where $g$ is the interaction strength, 
%$\sigma_z, \sigma_x$ are the Pauli matrices describing the two atomic levels separated 
%by energy difference $2\Delta$, and $b^\dagger$ ($b$) are 
%creation (annihilation) operators of a boson mode with frequency $\omega$. 
Let us make a canonical Bogoliubov transformation from $b, b^\dagger$ to the squeezed bosons $a, a^\dagger$ \cite{Emary02},
\beq
b=\frac{a+\tau\, a^\dagger}{\sqrt{1-\tau^2}}, ~~~~~~b^\dagger=\frac{\tau\,a+a^\dagger}{\sqrt{1-\tau^2}},
\eeq
where $|\tau|<1$ is a real parameter. In terms of the squeezed bosons, the Hamiltonian (\ref{2-photon-RabiH})
takes the form
\beqa
\tilde{H}_{2p}&=&\Delta\sigma_z+\frac{1}{1-\tau^2}\left[\left(\omega\tau+g\sigma_x(1+\tau^2)\right)\left((a^\dagger)^2 +a^2\right)\right.\n
 & &\left.+\left(\omega(1+\tau^2)+4g\tau\sigma_x\right)a^\dagger a+\omega\tau^2+2g\tau\sigma_x\right].
 \label{2-photon-RabiH1}
\eeqa
Introduce the operators $K_\pm, K_0$ 
\beq
K_+=\frac{1}{2}(a^\dagger)^2,~~~~K_-=\frac{1}{2}a^2,~~~~K_0=\frac{1}{2}\left(a^\dagger a+\frac{1}{2}\right).
   \label{2-photon-Rabi-boson}
\eeq
Then (\ref{2-photon-RabiH1}) becomes
%\begin{widetext}
%\beq
%\tilde{H}_{2-p}=\left(
%\begin{array}{cc}
%2\frac{\omega\tau+g(1+\tau^2)}{1-\tau^2} (K_++K_-) & \Delta \\
%+2\frac{\omega(1+\tau^2)+4g\tau}{1-\tau^2}K_0-\frac{1}{2}\omega} & \\
%\Delta & 2\frac{\omega\tau-g(1+\tau^2)}{1-\tau^2}(K_++K_-) \\
% & +2\frac{\omega(1+\tau^2)-4g\tau}{1-\tau^2}K_0-\frac{1}{2}\omega}
%\end{array}
%\right).
%\eeq
%\end{widetext}
\beqa
\tilde{H}_{2p}&=&\Delta\sigma_z+\frac{1}{1-\tau^2}\left[2\left(\omega\tau+g\sigma_x(1+\tau^2)\right)(K_++K_-)\right.\n
& &\left.+2\left(\omega(1+\tau^2)+4g\tau\sigma_x\right)K_0\right]-\frac{1}{2}\omega.
\eeqa

The operators $K_\pm, K_0$ form the $su(1,1)$ Lie algebra. Its quadratic Casimir,
$C=K_+K_--K_0(K_0-1)$, takes the particular values $C=\frac{3}{16}$ in the representation 
(\ref{2-photon-Rabi-boson}).  This is the well-known infinite-dimensional unitary irreducible
representation ${\cal D}^+(q)$ of $su(1,1)$ with $q=\frac{1}{4},  \frac{3}{4}$.
Thus the Fock-Hilbert space decomposes into the direct sum of two subspaces ${\cal H}^q$ 
labeled by $q=1/4, 3/4$.

%The basis state in the subspace ${\cal H}^q$, denoted as $|q,n\ra,~n=0,1,2,\cdots$, has the form
%\beq
%|q,n\ra=\frac{(a^\dagger)^{2(n+q-\frac{1}{4})}}{\sqrt{\left[2(n+q-\frac{1}{4}\right]!}}|0\ra,
%\eeq
%and the action of $\{K_\pm, K_0\}$ in this representation is given by
%\beqa
%K_0|q,n\ra&=&(n+q)|q,n\ra,\n
%K_+|q,n\ra&=&\sqrt{\left(n+q+\frac{1}{4}\right)\left(n+q+\frac{3}{4}\right)}|q,n+1\ra,\n
%K_-|q,n\ra&=&\sqrt{\left(n+q-\frac{1}{4}\right)\left(n+q-\frac{3}{4}\right)}|q,n-1\ra.\label{2-photon-rep}
%\eeqa
%Using the Fock-Bargmann correspondence  
%\beq
%a^{\dg} \longrightarrow z , \hspace{1cm} a \longrightarrow \frac{d}{dz},
%\hspace{1cm} |n\ra \longrightarrow \frac{z^{n}}{\sqrt{n!}},\label{Fock-Bargmann-correspondence}
%\eeq
%we can make the following association
%\beqa
%|q,n\ra  \longrightarrow \Psi_{q,n}(z)=\frac{z^n}{\sqrt{\left( 2(n+q-\frac{1}{4})\right)!}}.
%\eeqa
%Thus in each subspace ${\cal H}^q$ of the Bargmann space with basis vectors $\Psi_{q,n}(z)$, 
Using the representation $a^\dagger=w$ and $a=\frac{d}{dw}$, we can show
\cite{Zhang13b} that in the Bargmann space with basis vectors given by the monomials in $z=w^2$,
$\left\{{z^n}/{\sqrt{[2(n+q-1/k^2)]!}}\right\}$, 
the operators $K_\pm, K_0$ (\ref{2-photon-Rabi-boson}) have the  
single-variable 2nd order differential realization 
\beq
{K}_0 = z\frac{d}{dz}+ q, ~~~~{K}_+ = \frac{z}{2}, ~~~~~~
{K}_- =2z\frac{d^2}{dz^2}+4q\frac{d}{dz}.\label{2-photon-Rabi-diff}
\eeq
The Bargmann space is the Hilbert space of entire functions on the complex plane if the inner product
\beq
(f,g)=\int\,\overline{f(w)}\,g(w)\,d\mu(w)
\eeq
is finite for an appropriate measure $d\mu(w)$. Take
$d\mu(w)=\frac{1}{\pi}|w|^{4(q-1/4)}\,e^{-|w|^2}\,dx\,dy$. Then we can show 
\beq
(z^m, z^n)=(w^{2m}, w^{2n})=\left[2\left(n+q-{1}/{4}\right)\right]!\,\delta_{m n}.
\eeq
Thus the monomials $\left\{{z^n}/{\sqrt{[2(n+q-1/k^2)]!}}\right\}$ form an orthonormal basis 
of the Bargamann-Hilbert space. It is now not difficult to see that if $f(z)=\sum_{n=0}^\infty\,
c_n z^n$ then
\beq
||f||^2=\sum_{n=0}^\infty \,|c_n|^2\,\left[2\left(n+q-{1}/{4}\right)\right]!\label{entireness-test1}
\eeq
and $f(z)$ is entire if the sum on the right hand side converges.

Using this differential realization, working in a representation defined by $\sigma_x$ diagonal
%the 2-photon Rabi matrix Hamiltonian becomes
%\begin{widetext}
%\beq
%\tilde{H}_{2p}=\left(
%\begin{array}{cc}
%\frac{\omega\tau+g(1+\tau^2)}{1-\tau^2} \left(z+4z\frac{d^2}{dz^2}+8q\frac{d}{dz}\right) & \Delta \\
%+2\frac{\omega(1+\tau^2)+4g\tau}{1-\tau^2}\left(z\frac{d}{dz}+ q\right)-\frac{1}{2}\omega & \\
%\Delta & \frac{\omega\tau-g(1+\tau^2)}{1-\tau^2}\left(z+4z\frac{d^2}{dz^2}+8q\frac{d}{dz}\right) \\
% & +2\frac{\omega(1+\tau^2)-4g\tau}{1-\tau^2}\left(z\frac{d}{dz}+ q\right)-\frac{1}{2}\omega
%\end{array}
%\right). \label{2-photon-RabiH2}
%\eeq
%\end{widetext}
and choosing $\tau$ such that $\omega\tau+g(1+\tau^2)=0$, i.e.\footnote{The quadratic equation has two
roots $\tau=-\frac{\omega}{2g}(1\mp\Omega)$. The root given in (\ref{tau}) is real and obeys 
 $|\tau|<1$, the requirements of $\tau$ from the Bogoliubov transformation, provided
that $\left|\frac{2g}{\omega}\right|<1$.  This is seen as follows.
Assume $|\tau|=\left|\frac{\omega}{2g}\right|\left(1-\sqrt{1-\left(\frac{2g}{\omega}\right)^2}\right)
\geq 1$, i.e. $1-\left|\frac{2g}{\omega}\right|\geq \sqrt{1-\left(\frac{2g}{\omega}\right)^2}$. 
Then we would have 
$\sqrt{1-\left|\frac{2g}{\omega}\right|}\geq \sqrt{1+\left|\frac{2g}{\omega}\right|}$, 
which is impossible. }
\beq
\tau=-\frac{\omega}{2g}(1-\Omega),~~\Omega=\sqrt{1-\frac{4g^2}{\omega^2}},\label{tau}
\eeq
where $\left|\frac{2g}{\omega}\right|<1$, then the transformed 2-photon Rabi Hamiltonian becomes 
the matrix differential operator
% \begin{widetext}
\beq
\tilde{H}_{2p}=\left(
\begin{array}{cc}
2\omega\Omega\left(z\frac{d}{dz}+ q\right)-\frac{1}{2}\omega & \Delta  \\
\Delta & -\frac{8g}{\Omega}z\frac{d^2}{dz^2}+\frac{2}{\Omega}\left[\omega(2-\Omega^2)z-8gq\right]\frac{d}{dz} \\
 & -\frac{2g}{\Omega}z+\frac{2\omega(2-\Omega^2)q}{\Omega}-\frac{1}{2}\omega
\end{array}
\right). \label{2-photon-RabiH3}
\eeq
%\end{widetext}

In terms of two-component wavefuntion $\psi(z)=\left(\psi_+(z), \psi_-(z)\right)^T$,
the time-independent Schr\"odinger equation, $\tilde{H}_{2p}\psi(z)=E\psi(z)$,
yields a system of coupled differential equations,
%\begin{widetext}
\beqa
&&\left[2\omega\Omega\left(z\frac{d}{dz}+ q\right)-\frac{1}{2}\omega-E\right]\psi_++\Delta\psi_-=0,
     \label{2-photon-diff1}\\
&&\left[8gz\frac{d^2}{dz^2}+\left(-2\omega(2-\Omega^2)z+16gq\right)\frac{d}{dz}   \right.\n
&&~~~~~~~~~\left. 
+2gz-2\omega(2-\Omega^2)q+(\frac{1}{2}\omega+E)\Omega\right]\psi_-  -\Omega\Delta\,\psi_+=0. 
     \label{2-photon-diff2}
\eeqa
%\end{widetext}
This is a system of differential equations of Fuchsian type. 
Solutions to these equations must be analytic in the whole complex plane if
$E$ belongs to the spectrum of $\tilde{H}_{2p}$. So we are seeking solutions of the form
\beqa
\psi_+(z)=\sum_{n=0}^\infty\, {\cal K}_n^+(E)\, z^n, ~~~~
\psi_-(z)=\sum_{n=0}^\infty\, {\cal K}_n^-(E)\, z^n, \label{2-photon-series-solution}
\eeqa
which converge in the entire complex plane, i.e. solutions which are entire. 

Substituting (\ref{2-photon-series-solution}) into (\ref{2-photon-diff1}), we obtain
\beq
{\cal K}^+_n=\frac{\Delta}{E+\frac{1}{2}\omega-(2n+2q)\omega\Omega}{\cal K}_n^-.\label{2-photon-K+K-relation}
\eeq
%It follows,
%\beq
%\psi_+(z)=\sum_{n=0}^\infty\frac{\Delta}{E+\omega-(2n+2q)\omega\Omega}{\cal K}_n^-(E)\,z^n.
%\eeq
So ${\cal K}^+_n$ is not analytic in $E$ but has simple poles at
\beq
E=-\frac{1}{2}\omega+(2n+2q)\omega\Omega,~~~~n=0, 1, \cdots.\label{2-photon-exceptional}
\eeq
The energies (\ref{2-photon-exceptional}) appear for special values of model parameters 
\cite{Emary02,Zhang13b} and 
correspond to the exceptional solutions of the 2-photon Rabi model. If (\ref{2-photon-exceptional})
is satisfied, the infinite series expansions
(\ref{2-photon-series-solution}) truncate and reduce to polynomials in $z$ but only if the system
parameters satisfy certain constraints \cite{Zhang13b}. 
Majority part of the spectrum of the 2-photon Rabi model is regular 
for which (\ref{2-photon-exceptional}) is not satisfied. 
The regular spectrum of the model is given by the zeros of the transcendental function $F(E)$ obtained below.
Thus similar to the Rabi case, the spectrum of the 2-photon Rabi model consists of two parts, the regular and
the exceptional spectrum.

{}From (\ref{2-photon-diff2}), we obtain the 3-step recurrence relation for ${\cal K}_n^-$,
\beqa
&&{\cal K}_1^-+A_0\,{\cal K}_0^-=0,\n
&&{\cal K}_{n+1}^-+A_n\,{\cal K}_n^-+B_n\,{\cal K}_{n-1}^-=0,~~~n\geq 1,   \label{2-photon-3-step}
\eeqa
where
\beqa
A_n&=&\frac{1}{8g(n+1)(n+2q)}\left[-(2n+2q)\omega(2-\Omega^2)\right.\n
& &     \left.+\left(E+\frac{1}{2}\omega-\frac{\Delta^2}{E+\frac{1}{2}\omega-(2n+2q)\omega\Omega}\right)\Omega\right],\n
B_n&=&\frac{1}{4(n+1)(n+2q)}.
\eeqa
The coefficients $A_n, B_n$ have the behavior as $n\rightarrow\infty$
\beq
A_n\sim a\,n^\alpha,~~~~~~B_n\sim b\,n^\beta
\eeq
with 
\beq
a=-\frac{\omega}{4g}(2-\Omega^2),~~~~\alpha=-1,~~~~b=\frac{1}{4},~~~~\beta=-2.\label{ab-values}
\eeq
Thus the asymptotic structure of solutions to the 2nd equation of (\ref{2-photon-3-step}) 
depends on the Newton-Puiseux
diagram formed with the points $P_0(0,0), P_1(1,-1), P_2(2,-2)$ \cite{Gautschi67}.
Let $\gamma$ be the slope of $\overline{P_0P_1}$ and $\delta$ the slope of $\overline{P_1P_2}$
so that $\gamma=\alpha$ and $\delta=\beta-\alpha$. Then we have $\gamma=\delta=\alpha$. The
characteristic equation of the $n\geq 1$ part of (\ref{2-photon-3-step}) reads 
$t^2+at+b=0$ with $a, b$ given in (\ref{ab-values}). It has two roots $t_1=\frac{\omega}{4g}$,
$t_2=\frac{g}{\omega}$. Remembering the condition $\left|\frac{2g}{\omega}\right|<1$, we have $|t_2|<|t_1|$.
Applying the Perron-Kreuser theorem (i.e. Theorem 2.3 of \cite{Gautschi67}), 
we conclude that the two linearly independent solutions ${\cal K}_{n,1}^-$ and ${\cal K}_{n,2}^-$
of the $n\geq 1$ part (i.e. the truly 3-term part) of (\ref{2-photon-3-step}) satisfy
\beq
\lim_{n\rightarrow\infty}\frac{{\cal K}_{n+1,r}^-}{{\cal K}_{n,r}^-}\sim t_r\,n^{-1},~~~~r=1,2.
   \label{K-asymptotics}
\eeq
So ${\cal K}_{n,2}^-$ is a minimal solution and ${\cal K}_{n,1}^-$ is a dominant one.
By (\ref{entireness-test1}), we can see that the infinite power series in 
(\ref{2-photon-series-solution}) with expansion coefficients ${\cal K}_{n,r}^-$ is entire if the sum
\beq
\sum_{n=0}^\infty \left|{\cal K}_{n,r}^-\right|^2[2(n+q-1/4)]!\label{K-sum}
\eeq
converges. Using the asymptotic form (\ref{K-asymptotics}) we get
\beq
\lim_{n\rightarrow\infty}\frac{\left|{\cal K}_{n+1,r}^-\right|^2[2(n+1+q-1/4)]!}
{\left|{\cal K}_{n,r}^-\right|^2[2(n+q-1/4)]!}=4\,|t_r|^2
\eeq
which is less than 1 for $r=2$ and greater than 1 for $r=1$. Thus by the ratio test, the sum
(\ref{K-sum}) converges for the minimal solution ${\cal K}_n^{min}\equiv {\cal K}_{n,2}^-$
and diverges for the dominant solution ${\cal K}_{n,1}^-$. It follows that
the infinite power series expansions $\psi_\pm^{min}(z)$, obtained
by substituting ${\cal K}_n^{min}$ for the ${\cal K}^-_n$'s in (\ref{2-photon-K+K-relation}) and 
(\ref{2-photon-series-solution}), converge in the whole complex plane, i.e. they are entire.

We now find energy eigenvalues $E$ corresponding to the
the minimal solution ${\cal K}_n^{min}$ (and thus to the entire wavefunctions $\psi_\pm^{min}(z)$).
We follow a procedure presented in \cite{Leaver86} that uses the relationship between minimal
solutions and infinite continued fractions \cite{Gautschi67}. Such a procedure was applied in 
\cite{Moroz12} to analyze the Rabi model.

By the Pincherle theorem (i.e. Theorem 1.1 of \cite{Gautschi67}), the ratio of successive elements of
the minimal solution sequence ${\cal K}_n^{min}$ is expressible as continued fractions,
%Proceeding in the direction of increasing $n$, we have
\beq
R_{n}=\frac{{\cal K}_{n+1}^{min}}{{\cal K}_n^{min}}=-\frac{B_{n+1}}{~A_{n+1}-}\,\frac{B_{n+2}}{~A_{n+2}-}\,
\frac{B_{n+3}}{~A_{n+3}-}\,\cdots,   \label{2-photon-continued-fraction}
\eeq
which for $n=0$ gives
\beq
R_{0}=\frac{{\cal K}_{1}^{min}}{{\cal K}_0^{min}}=-\frac{B_{1}}{~A_{1}-}\,\frac{B_{2}}{~A_{2}-}\,
\frac{B_{3}}{~A_{3}-}\,\cdots.   \label{2-photon-continued-fraction1}
\eeq
Note that the ratio $R_0=\frac{{\cal K}_1^{min}}{{\cal K}_0^{min}}$ involves ${\cal K}_0^{min}$, 
although the above continued fraction expression is obtained from the 
2nd equation of (\ref{2-photon-3-step}), i.e the recurrence (\ref{2-photon-3-step}) for $n\geq 1$. 
However, for single-ended sequences such as those appearing in the infinite series expansions
(\ref{2-photon-series-solution}), the ratio $R_0=\frac{{\cal K}_1^{min}}{{\cal K}_0^{min}}$ of 
the first two terms of a minimal solution is unambiguously fixed by the first equation of 
the recurrence (\ref{2-photon-3-step}), namely,
%\begin{widetext}
\beq
 R_0=-A_0=\frac{1}{16gq}\left[2q\omega(2-\Omega^2)
   -\left(E+\frac{1}{2}\omega-\frac{\Delta^2}{E+\frac{1}{2}\omega-2q\omega\Omega}\right)\Omega\right]   
   \label{2-photon-continued-fraction2}
\eeq
%\end{widetext}
In general, the $R_0$ computed from the continued fraction (\ref{2-photon-continued-fraction1}) 
can not be the same as that from (\ref{2-photon-continued-fraction2})
%i.e. (\ref{2-photon-continued-fraction1}) and (\ref{2-photon-continued-fraction2}) are not both satisfied, 
for arbitrary values of
recurrence coefficients $A_n$ and $B_n$. As a result, general solutions to the recurrence 
(\ref{2-photon-3-step}) are dominant and are usually
generated by simple forward recursion from a given value of ${\cal K}_0^-$. Physical meaningful
solutions are those that are entire in the Bargmann-Hilbert space \cite{Schweber67}. 
They can be obtained if $E$ can be adjusted so that equations 
(\ref{2-photon-continued-fraction1}) and (\ref{2-photon-continued-fraction2})
are both satisfied. Then the resulting solution sequence ${\cal K}_n^-(E)$ will be purely minimal and
the power series expansion (\ref{2-photon-series-solution}) will converge in the whole complex plane.

Therefore, if we define the transcendental function $F(E)=R_0+A_0$ with $R_0$  
given by the continued fraction in (\ref{2-photon-continued-fraction1}), then the
 zeros of $F(E)$ correspond to the points in the parameter space where the condition 
(\ref{2-photon-continued-fraction2}) is satisfied. In other words, $F(E)=0$ is the eigenvalue equation
of the 2-photon Rabi model, which may be solved for $E$ by standard nonlinear root-search techniques
(see e.g. \cite{Leaver86,Liu92} and references therein).
 Only for the denumerable infinite values of $E$ which are the roots of $F(E)=0$, do
we get entire solutions of the differential equations (\ref{2-photon-diff1}) and (\ref{2-photon-diff2}).



\section{Two-mode quantum Rabi model}\label{two-mode}
The Hamiltonian of the two-mode Rabi model is
\beq
H_{2m}=\omega(b_1^\dagger b_1+b_2^\dagger b_2)+\Delta\sigma_z+g\,\sigma_x(b_1^\dagger b_2^\dagger +b_1 b_2),
  \label{2-mode-RabiH}
\eeq
where we assume that the boson modes are degenerate with the same frequency $\omega$.
Introduce the two-mode Bogoliubov transformation,
\beqa
b_1=\frac{a_1+\sigma\, a_2^\dagger}{\sqrt{1-\sigma^2}},~~~~
  b_1^\dagger=\frac{\sigma\, a_2+a_1^\dagger}{\sqrt{1-\sigma^2}},~~~~ %\n
b_2=\frac{a_2+\sigma\, a_1^\dagger}{\sqrt{1-\sigma^2}},~~~~
  b_2^\dagger=\frac{\sigma\, a_1+a_2^\dagger}{\sqrt{1-\sigma^2}}.
\eeqa
Here $|\sigma|<1$ is a real parameter and $a_1, a_2, a_1^\dagger, a_2^\dagger$
are squeezed bosons satisfying the canonical commutation relations 
%\beq
$[a_i, a_i^\dagger]=1,~~[a_i, a_j]=[a_i, a_j^\dagger]=[a_i^\dagger, a_j^\dagger]=0,~i,j=1,2$.
%\eeq
In terms of the 2-mode squeezed bosons, the Hamiltonian (\ref{2-mode-RabiH}) has the form
%\begin{widetext}
\beqa
\tilde{H}_{2m}&=&\frac{1}{{1-\sigma^2}}\left[\left(2\omega\sigma+g\sigma_x(1+\sigma^2)\right)
    (a_1^\dagger a_2^\dagger+a_1 a_2)\right.\n
& &+\left(\omega(1+\sigma^2)+2g\sigma\,\sigma_x\right)(a_1^\dagger a_1+a_2^\dagger a_2) %\n
\left. +2\omega\sigma^2+2g\sigma\,\sigma_x\right]+\Delta\sigma_z. \label{2-mode-RabiH1}
\eeqa
%\end{widetext}
Introduce the operators $K_\pm, K_0$
\beq
K_+=a_1^\dagger a_2^\dagger,~~~~K_-=a_1 a_2,~~~~K_0=\frac{1}{2}(a_1^\dagger a_1+a_2^\dagger a_2+1).\label{2-mode-boson}
\eeq
Then (\ref{2-mode-RabiH1}) becomes
%\begin{widetext}
%\beq
%\tilde{H}_{2m}=\left(
%\begin{array}{cc}
%\frac{2\omega\sigma+g(1+\sigma^2)}{1-\sigma^2} (K_++K_-) & \Delta \\
%+2\frac{\omega(1+\sigma^2)+2g\sigma}{1-\sigma^2}K_0-\omega & \\
%\Delta & \frac{2\omega\sigma-g(1+\sigma^2)}{1-\sigma^2} (K_++K_-) \\
% & +2\frac{\omega(1+\sigma^2)-2g\sigma}{1-\sigma^2}K_0-\omega
%\end{array}
%\right). \label{2-mode-RabiH2}
%\eeq
%\end{widetext}
\beqa
\tilde{H}_{2m}&=&\Delta\sigma_z+\frac{1}{{1-\sigma^2}}\left[\left(2\omega\sigma+g\sigma_x(1+\sigma^2)\right)
    (K_++K_-)\right.\n
& &\left. +2\left(\omega(1+\sigma^2)+2g\sigma\sigma_x\right)K_0\right] -\omega.\label{2-mode-RabiH2}
\eeqa

The operators $K_\pm, K_0$ form the $su(1,1)$ Lie algebra. Its quadratic Casimir, $C=K_+K_--K_0(K_0-1)$, 
takes the particular values $C=\kappa(1-\kappa)$ in the representation (\ref{2-mode-boson}), 
where $\kappa=1/2, 1, 3/2,\cdots$. This is the well-known infinite-dimensional unitary 
irreducible representation of $su(1,1)$ known as the positive discrete series 
${\cal D}^+(\kappa)$. Thus the Fock-Hilbert space decomposes 
into the direct sum of infinite subspaces ${\cal H}^\kappa$ labeled by $\kappa=1/2, 1, 3/2, \cdots$.

%The basis state in the subspace ${\cal H}^\kappa$, denoted as $|\kappa,n\ra,~n=0,1,2,\cdots$, 
%has the form
%\beq
%|\kappa,n\ra=\frac{(a_1^\dagger)^{n+2\kappa-1} (a_2^\dagger)^n}{\sqrt{(n+2\kappa-1)!n!}}|0\ra,
%\eeq
%and the action of $K_\pm, K_0$ in this representation is given by
%\beqa
%K_0|\kappa,n\ra&=&(n+\kappa)|\kappa,n\ra,\n
%K_+|\kappa,n\ra&=&\sqrt{(n+2\kappa)(n+1)}|\kappa,n+1\ra,\n
%K_-|\kappa,n\ra&=&\sqrt{(n+2\kappa-1)n}|\kappa,n-1\ra.\label{2-mode-rep}
%\eeqa

%Using the Fock-Bargmann correspondence (\ref{Fock-Bargmann-correspondence}),
%we can show that the infinite set of monomials
%\beq
%\Psi_{\kappa,n}(z)=\frac{z^n}{\sqrt{(n+2\kappa-1)!n!}},~~~~n=0,1,2,\cdots,
%\eeq
%form the basis in the Bargmann-Hilbert subspace associated with the representation (\ref{2-mode-rep}). 
Using the representation $a_i^\dagger=w_i$ and $a_i=\frac{d}{dw_i}$, $i=1,2$, we may show 
\cite{Zhang13b} that in the Bargmann space with basis vectors given by the monomials in $z=w_1 w_2$,
$\left\{{z^n}/{\sqrt{n!(n+2\kappa-1)!}}\right\}$, the operators $K_\pm, K_0$ (\ref{2-mode-boson}) 
have single-variable differential realization
\beq
K_0=z\frac{d}{dz}+\kappa,~~~~K_+=z,~~~~K_-=z\frac{d^2}{dz^2}+2\kappa\frac{d}{dz},
      \label{su11-diff-rep-2mode}
\eeq
where $\kappa=1/2, 1, 3/2,\cdots$.
The Bargmann space is the Hilbert space of entire functions on ${\bf C}^2$ if the inner product
\beq
(f,g)=\int\int\,\overline{f(w_1,w_2)}\,g(w_1,w_2)\,d\mu(w_1,w_2)
\eeq
is finite for an appropriate measure $d\mu(w_1,w_2)$. It can be easily seen that if we choose
$d\mu(w_1, w_2)=\frac{1}{\pi^2}|w_1|^{2(2\kappa-1)}\,e^{-|w_1|^2-|w_2|^2}\,d^2w_1\,d^2w_2$,
then 
\beq
(z^m, z^n)=((w_1w_2)^{2m}, (w_1w_2)^{2n})=n!\left(n+2\kappa-1\right)!\,\delta_{m n}.
\eeq
Thus the monomials $\left\{{z^n}/{\sqrt{n!(n+2\kappa-1)!}}\right\}$ form an orthonormal basis 
of the Bargamann-Hilbert space. It is now not difficult to see that if $f(z)=\sum_{n=0}^\infty\,
c_n z^n$ then
\beq
||f||^2=\sum_{n=0}^\infty \,|c_n|^2\,n!\left(n+2\kappa-1\right)!\label{entireness-test2}
\eeq
and $f(z)$ is entire if the sum on the right hand side converges.

Using this differential realization, working in a representation defined by $\sigma_x$ diagonal
% we can express the 
%Hamiltonian (\ref{2-mode-RabiH2}) as the 2nd-order matrix differential operator
%in each Bargmann-Hilbert subspace labeled by $\kappa$,
%\begin{widetext}
%\beq
%\tilde{H}_{2m}=\left(
%\begin{array}{cc}
%\frac{2\omega\sigma+g(1+\sigma^2)}{1-\sigma^2} \left(z+z\frac{d^2}{dz^2}+2\kappa\frac{d}{dz}\right) & \Delta \\
%+2\frac{\omega(1+\sigma^2)+2g\sigma}{1-\sigma^2}\left(z\frac{d}{dz}+\kappa\right)-\omega & \\
%\Delta & \frac{2\omega\sigma+g(1+\sigma^2)}{1-\sigma^2} \left(z+z\frac{d^2}{dz^2}+2\kappa\frac{d}{dz}\right)\\
% & +2\frac{\omega(1+\sigma^2)+2g\sigma}{1-\sigma^2}\left(z\frac{d}{dz}+\kappa\right)-\omega
%\end{array}
%\right).
%\eeq
%\end{widetext}
and choosing $\sigma$ such that $2\omega\sigma+g(1+\sigma^2)=0$, i.e.\footnote{By arguments similar to
the 2-photon case we can show that $\sigma$ given by (\ref{sigma}) obeys the requirements that 
it is real and $|\sigma|<1$, provided that $\left|\frac{g}{\omega}\right|<1$.}
\beq
\sigma=-\frac{\omega}{g}(1-\Lambda),~~~~\Lambda=\sqrt{1-\frac{g^2}{\omega^2}},\label{sigma}
\eeq
where $\left|\frac{g}{\omega}\right|<1$, then the transformed  Hamiltonian (\ref{2-mode-RabiH2})
 becomes a matrix differential operator
%\begin{widetext}
\beq
\tilde{H}_{2m}=\left(
\begin{array}{cc}
2\omega\Lambda\left(z\frac{d}{dz}+\kappa\right)-\omega & \Delta \\
\Delta & -\frac{2g}{\Lambda}z\frac{d^2}{dz^2}+\frac{2}{\Lambda}\left[\omega(2-\Lambda^2)z-2g\kappa\right]\frac{d}{dz} \\
 & -\frac{2g}{\Lambda}z+\frac{2\omega(2-\Lambda^2)\kappa}{\Lambda}-\omega
\end{array}
\right).\label{2-mode-RabiH3}
\eeq
%\end{widetext}
In terms of two-component wavefuntion $\phi(z)=\left(\phi_+(z), \phi_-(z)\right)^T$, 
the time-independent Schr\"odinger equation, $\tilde{H}_{2m}\phi(z)=E\phi(z)$,
yields a system of coupled differential equations,
%\begin{widetext}
\beqa
&&\left[2\omega\Lambda\left(z\frac{d}{dz}+\kappa\right)-\omega-E\right]\phi_++\Delta\phi_-=0,
   \label{2-mode-diff1}\\
&&\left[2gz\frac{d^2}{dz^2}+\left(-2\omega(2-\Lambda^2)z+4g\kappa\right)\frac{d}{dz}  \right.\n
&&~~~~~~~~~\left.
+2gz-2\omega(2-\Lambda^2)\kappa+(E+\omega)\Lambda\right]\phi_- -\Lambda\Delta\,\phi_+=0.
   \label{2-mode-diff2}
\eeqa
%\end{widetext}
This is a system of differential equations of Fuchsian type. 
Solutions to these equations must be analytic in the whole complex plane if
$E$ belongs to the spectrum of $\tilde{H}_{2m}$. Similar to the 2-photon Rabi case, we seek solutions of the form
\beqa
\phi_+(z)=\sum_{n=0}^\infty\, {\cal Q}_n^+(E)\, z^n, ~~~~~~
\phi_-(z)=\sum_{n=0}^\infty\, {\cal Q}_n^-(E)\, z^n, \label{2-mode-series-solution}
\eeqa
which converge in the entire complex plane.
% i.e. solutions $\phi_\pm(z)$ which are entire in the Bargmann-Hilbert space. 

Substituting (\ref{2-mode-series-solution}) into (\ref{2-mode-diff1}), we obtain
\beq
{\cal Q}^+_n=\frac{\Delta}{E+\omega-(2n+2\kappa)\omega\Lambda}{\cal Q}_n^-.\label{2-mode-K+K-relation}
\eeq
%It follows,
%\beq
%\phi_+(z)=\sum_{n=0}^\infty\frac{\Delta}{E+\omega-(2n+2\kappa)\omega\Lambda}{\cal Q}_n^-(E)\,z^n.
%\eeq
So ${\cal Q}^+_n$ is not analytic in $E$ but has simple poles at
\beq
E=-\omega+(2n+2\kappa)\omega\Lambda,~~~~n=0, 1, \cdots.\label{2-mode-exceptional}
\eeq
The energies (\ref{2-mode-exceptional}) appear for special values of model parameters 
\cite{Zhang13b} and correspond to the
exceptional solutions of the two-mode Rabi model. If (\ref{2-mode-exceptional}) is satisfied,
the infinite series expansions (\ref{2-mode-series-solution}) truncate 
and reduce to polynomials in $z$ but only if the model parameters obey certain constraints \cite{Zhang13a}.
Majority part of the spectrum of the two-mode Rabi model is regular spectrum which does not have the form 
(\ref{2-mode-exceptional}). The regular spectrum of the model is given by the zeros of 
the transcendental function $G(E)$ obtained below. 
Thus again, the spectrum of the two-mode Rabi model consists of
two parts, the regular and the exceptional spectrum.

{}From (\ref{2-mode-diff2}), we obtain the 3-step recurrence relation for ${\cal Q}_n^-$,
\beqa
&&{\cal Q}_1^-+C_0\,{\cal Q}_0^-=0,\n
&&{\cal Q}_{n+1}^-+C_n\,{\cal Q}_n^-+D_n\,{\cal Q}_{n-1}^-=0,~~~n\geq 1, \label{2-mode-3-step}
\eeqa
where
\beqa
C_n&=&\frac{1}{2g(n+1)(n+2\kappa)}\left[-(2n+2\kappa)\omega(2-\Lambda^2)\right.\n
& &     \left.+\left(E+\omega-\frac{\Delta^2}{E+\omega-(2n+2\kappa)\omega\Lambda}\right)\Lambda\right],\n
D_n&=&\frac{1}{(n+1)(n+2\kappa)}.
\eeqa
The coefficients $C_n, D_n$ have the behavior as $n\rightarrow\infty$
\beq
C_n\sim c\,n^\mu,~~~~~~D_n\sim d\,n^\rho
\eeq
with 
\beq
c=-\frac{\omega}{g}(2-\Lambda^2),~~~~\mu=-1,~~~~d=1,~~~~\rho=-2.
\eeq
%The asymptotic structure of solutions to (\ref{2-mode-3-step}) depends on the Newton-Puiseux
%diagram formed with the points $P_0(0,0), P_1(1,-1), P_2(2,-2)$ \cite{Gautschi67}.
%Let $\mu$ be the slope of $\overline{P_0P_1}$ and $\nu$ the slope of $\overline{P_1P_2}$
%so that $\mu=\lambda$ and $\nu=\rho-\lambda$. Then we have $\mu=\nu=\lambda$. The
%characteristic equation of the $n\geq 1$ part of (\ref{2-mode-3-step}), 
%$t^2=0$, has two equal solutions $t_1=t_2=0$.
%Applying the Perron-Kreuser theorem (i.e. Theorem 2.3 of \cite{Gautschi67}),
By analysis similar to the 2-photon case,  we see that the two linearly independent solutions 
${\cal Q}_{n,1}^-$ and ${\cal Q}_{n,2}^-$ of the $n\geq 1$ part
 of the recurrence (\ref{2-photon-3-step}) obey
\beq
\lim_{n\rightarrow\infty}\frac{{\cal Q}_{n+1,r}^-}{{\cal Q}_{n,r}^-}\sim t_r\,n^{-1},~~~~r=1,2,
\eeq
where $t_1=\frac{\omega}{g}, ~ t_2=\frac{g}{\omega}$ and $|t_2|<|t_1|$ (from the condition 
$\left|\frac{g}{\omega}\right|<1$). Thus ${\cal Q}_{n,2}^-$ is a minimal solution and 
${\cal Q}_{n,1}^-$ is a dominant one. Using (\ref{entireness-test2}) and by similar analysis to
the 2-photon case, we can conclude that 
the infinite power series expansions $\phi^{min}_\pm(z)$
generated by substituting the minimal solution ${\cal Q}_n^{min}\equiv {\cal Q}_{n,2}^-$ 
for the ${\cal Q}_n^-$'s in 
(\ref{2-mode-K+K-relation}) and (\ref{2-mode-series-solution}), converge in the whole complex plane. 

%By the Pincherle theorem (Theorem 1.1 of \cite{Gautschi67}), 
The ratio of successive elements
of the minimal solution ${\cal Q}_n^{min}$ can be expressed as continued fractions,
%Proceeding in the direction of increasing $n$, we find
\beq
S_{n}=\frac{{\cal Q}_{n+1}^{min}}{{\cal Q}_n^{min}}=-\frac{D_{n+1}}{~C_{n+1}-}\,\frac{D_{n+2}}{~C_{n+2}-}\,
\frac{D_{n+3}}{~C_{n+3}-}\,\cdots,   \label{2-mode-continued-fraction}
\eeq
which for $n=0$ reduces to
\beq
S_{0}=\frac{{\cal Q}_{1}^{min}}{{\cal Q}_0^{min}}=-\frac{D_{1}}{~C_{1}-}\,\frac{D_{2}}{~C_{2}-}\,
\frac{D_{3}}{~C_{3}-}\,\cdots.   \label{2-mode-continued-fraction1}
\eeq
The ratio $S_0=\frac{{\cal Q}_1^{min}}{{\cal Q}_0^{min}}$ involves ${\cal Q}_0^{min}$, 
although the above continued fraction expression is obtained from the 2nd equation 
of (\ref{2-mode-3-step}). % i.e the recurrence (\ref{2-mode-3-step}) for $n\geq 1$. 
On the other hand, 
%for single-ended sequences such as those appearing in the infinite series expansions (\ref{2-mode-series-solution}), 
the ratio $S_0=\frac{{\cal Q}_1^{min}}{{\cal Q}_0^{min}}$ of the first two terms of 
a minimal solution is unambiguously fixed by  the $n=0$ part of 
the recurrence (\ref{2-mode-3-step}), that is,
%\begin{widetext}
\beq
S_0=-C_0=\frac{1}{4g\kappa}\left[2\kappa\omega(2-\Lambda^2)
   -\left(E+\omega-\frac{\Delta^2}{E+\omega-2\kappa\omega\Lambda}\right)\Lambda\right]   
   \label{2-mode-continued-fraction2}
\eeq
%\end{widetext}
In general, (\ref{2-mode-continued-fraction1})  and (\ref{2-mode-continued-fraction2}) can not be both satisfied
% i.e. the $S_0$ computed from (\ref{2-mode-continued-fraction1}) is not the same as that from 
%(\ref{2-mode-continued-fraction2}) 
for arbitrary values of the recurrence coefficients $C_n$ and $D_n$.  
Thus general solutions to the recurrence (\ref{2-mode-3-step}) are dominant and are usually
generated by simple forward recursion from a given value of ${\cal Q}_0^-$. 
Physical meaningful solutions are those that are entire in the Bargmann-Hilbert space
\cite{Schweber67}. They can be obtained if $E$ can be adjusted so that equations 
(\ref{2-mode-continued-fraction1}) and (\ref{2-mode-continued-fraction2}) are both satisfied. 
Then the resulting solution sequence ${\cal Q}_n^-(E)$ will be purely minimal and
the power series expansion (\ref{2-mode-series-solution}) will converge in the whole complex plane.
Equating the right hand sides of (\ref{2-mode-continued-fraction1}) and (\ref{2-mode-continued-fraction2})
yields an implicit continued fraction equation for the regular spectrum $E$. That is,
the regular energies $E$ of the two-mode Rabi model are determined by the zeros of the transcendental 
function $G(E)=S_0+C_0$ with $S_0$ and $C_0$ given by
(\ref{2-mode-continued-fraction1}) and (\ref{2-mode-continued-fraction2}), respectively.
Only for the denumerable infinite values of $E$ which are the roots of $G(E)=0$, do
we get entire solutions of the differential equations (\ref{2-mode-diff1}) and (\ref{2-mode-diff2}).
The transcendental eigenvalue equation $G(E)=0$ may be solved for $E$ by standard root-search algorithms
(see e.g. \cite{Leaver86,Liu92} and references therein). 


\section{Conclusions}\label{summary}
We have presented analytic representations of solutions to the 2-photon and two-mode quantum Rabi models.
Our main result is that the respective spectrum of the models 
consists of two parts, the regular and exceptional spectrum. Exceptional solutions of the models
appear for special values of system parameters and were previously found in \cite{Emary02,Zhang13b}.  
Majority of the eigenvalues of the models are regular and correspond to arbitrary coupling strength 
and detuning. The regular eigenvalues of the models are given by the zeros of the transcendental 
functions, which have been analytically found by applying the Bogoliubov
transformations and Bargmann-Hilbert spaces.  The zeros can be determined numerically by standard 
root-search techniques. 

This work should be of general interest due to its focus on the two simplest generalizations of
the popular quantum Rabi model, but also be intriguing to those interested in analytic solutions 
of Fuchsian differential equations in physics.
Our approach is quite general and could be applied to a large class of dynamical systems. As an example,
in the Appendix, we will solve the driven Rabi model with broken ${\cal Z}_2$ symmetry \cite{Braak11}.




%\vskip.3in
%\noindent{\bf Acknowledgments:}
\section*{Acknowledgments}
We would like to thank Alexander Moroz for pointing out an error in the previous version and for 
some useful comments and discussions. This work was partially supported by the Australian Research Council 
through Discovery-Projects grants DP110103434 and DP140101492. 





\appendix
\section{Driven Rabi model}
In this appendix we revisit the driven quantum Rabi model
and derive a transcendental equation for its regular energy spectrum.


The driven Rabi model is a generalization of the Rabi model
with broken ${\cal Z}_2$ symmetry. It has recently been used in \cite{Larson13} to examine
quantum thermalization. The Hamiltonian of the model reads
\beq
H_{dR}=\omega b^\dagger b+\Delta\,\sigma_z+g\,\sigma_x(b^\dagger+b)+\delta\,\sigma_x,\label{RabiH1}
\eeq
where $\delta$ is the drive amplitude. 
After a canonical Bogoliubov transformation $b=a+\frac{\lambda}{\omega},~~b^\dagger=a^\dagger+\frac{\lambda}{\omega}$, 
where $\lambda$ is a real parameter, the Hamiltonian (\ref{RabiH1}) becomes
\beq
\tilde{H}_{dR}=\omega a^\dagger a+\Delta\,\sigma_z+(g\,\sigma_x+\lambda) (a^\dagger+a)
   +\delta\,\sigma_x+\frac{2g\lambda}{\omega}\,\sigma_x+\frac{\lambda^2}{\omega}.
\eeq

Using the Bargmann realization $a^\dagger=z,~~a=\frac{d}{dz}$, 
working in a representation defined by $\sigma_x$ diagonal and choosing $\lambda=-g$,
we can turn the transformed driven Rabi Hamiltonian into the matrix differential operator
\beq
\tilde{H}_{dR}=\left(
\begin{array}{cc}
\omega\,z\frac{d}{dz}+\delta-\frac{g^2}{\omega} & \Delta \\
\Delta & \omega\left(z-\frac{2g}{\omega}\right)\frac{d}{dz}-2gz-\delta +\frac{3g^2}{\omega}
\end{array}
\right ).
\eeq
In this case, the Bargmann-Hilbert space is the space of entire functions $f(z)$ with inner product
(\ref{inner-product1}) and orthonormal basis $\left\{\frac{z^n}{\sqrt{n!}}\right\}$. So if
$f(z)=\sum_{n=0}^\infty c_n\,z^n$, then 
\beq
||f||^2=\sum_{n=0}^\infty|c_n|^2\,n!
\eeq
and $f(z)$ is entire iff this sum converges.

In terms of two-component wavefunction $\varphi(z)=(\varphi_+(z), \varphi_-(z))^T$, the
time-independent Schr\"odinger equation, $\tilde{H}_{dR}\,\varphi(z)=E\,\varphi(z)$, yields a system
of couple of differential equations
\beqa
&&\left(\omega\,z\frac{d}{dz}+\delta-\frac{g^2}{\omega}-E\right)\varphi_+ +\Delta\,\varphi_- =0,\label{Rabi-diff1}\\
&&\left[\omega\left(z-\frac{2g}{\omega}\right)\frac{d}{dz}-2gz+\frac{3g^2}{\omega}-\delta-E\right]\varphi_-
    +\Delta\,\varphi_+ =0.\label{Rabi-diff2}
\eeqa
Solutions to these equations must be analytic in the whole complex plane 
if $E$ belongs to the spectrum of $\tilde{H}_R$. So we are looking for
solutions of the form
\beqa
\varphi_+(z)=\sum_{n=0}^\infty {\cal R}_n^+(E)\,z^n,~~~~~~
\varphi_-(z)=\sum_{n=0}^\infty {\cal R}_n^-(E)\,z^n,\label{Rabi-series-solution}
\eeqa
which converge in the entire complex plane, i.e. solutions which are entire.

Substituting (\ref{Rabi-series-solution}) into (\ref{Rabi-diff1}), we obtain
\beq
{\cal R}_n^+=\frac{\Delta}{E-n\omega-\delta+\frac{g^2}{\omega}}\,{\cal R}_n^-.
\eeq
Thus ${\cal R}_n^+$ is not analytic in $E$ but has simple poles at
\beq
E=n\omega+\delta-\frac{g^2}{\omega},~~~~n=0,1,\cdots.
   \label{Rabi-exceptional}
\eeq
The energies (\ref{Rabi-exceptional}) appear for special value of model parameters and 
correspond to the exceptional solutions of the driven Rabi model. 
If (\ref{Rabi-exceptional}) is satisfied, the infinite series
expansions (\ref{Rabi-series-solution}) truncate and reduce to polynomials in $z$ but only
if the system parameters satisfy certain constraints. This can be easily verified by following
the procedure in \cite{Zhang13b} (and thus the driven Rabi model is quasi-exactly solvable).  
Majority part of 
the spectrum of the driven Rabi model is regular for which (\ref{Rabi-exceptional}) is not satisfied.
The regular spectrum of the model is given by the zeros of the transcendental function 
$Q(E)$ obtained below.

{}From (\ref{Rabi-diff2}), we obtain the 3-term recurrence relation for ${\cal R}_n^-$
\beqa
&&{\cal R}_1^-+X_0\,{\cal R}_0^-=0,\n
&&{\cal R}_{n+1}^-+X_n\,{\cal R}_n^-+Y_n\,{\cal R}_{n-1}^-=0,~~~n\geq 1,\label{Rabi-3term}
\eeqa
where
\beqa
&&X_n=\frac{1}{2g(n+1)}\left[E-n\omega+\delta-\frac{3g^2}{\omega}+\frac{\Delta^2}
   {E-n\omega-\delta+\frac{g^2}{\omega}}\right],\n
&&Y_n=\frac{1}{n+1}.
\eeqa
The characteristic equation for the $n\geq 1$ part of the recurrence relation (\ref{Rabi-3term}) is 
$t^2-\frac{\omega}{2g}\,t=0$, which gives two distinct roots $t_1=0$ and $t_2=\frac{\omega}{2g}$.
Applying the Poincar\'e-Perron theorem (i.e. theorems 2.1 and 2.2 of \cite{Gautschi67}), the two
linearly independent solutions ${\cal R}_{n, 1}^-$ and ${\cal R}_{n, 2}^-$ of 
the 2nd equation of (\ref{Rabi-3term}) satisfy
\beq
\lim_{n\rightarrow\infty}\,\frac{{\cal R}_{n+1, r}^-}{{\cal R}_{n, r}^-}=t_r,~~~~~r=1, 2,
\eeq
and ${\cal R}_{n}^{min}\equiv {\cal R}_{n, 1}^-$ is a minimal solution. 
It is not difficult to show that the infinite series expansions (\ref{Rabi-series-solution}) with
coefficients given by the minimal solution ${\cal R}_{n}^{min}$ are entire.

The ratio of successive elements of the minimal solution sequence ${\cal R}_{n}^{min}$ 
is expressible in terms of continued fractions,
\beq
T_{n}=\frac{{\cal R}_{n+1}^{min}}{{\cal R}_n^{min}}=-\frac{Y_{n+1}}{~X_{n+1}-}\,\frac{Y_{n+2}}{~X_{n+2}-}\,
\frac{Y_{n+3}}{~X_{n+3}-}\,\cdots,   \label{Rabi-continued-fraction}
\eeq
which for $n=0$ reduces to
\beq
T_{0}=\frac{{\cal R}_{1}^{min}}{{\cal R}_0^{min}}=-\frac{Y_{1}}{~X_{1}-}\,\frac{Y_{2}}{~X_{2}-}\,
\frac{Y_{3}}{~X_{3}-}\,\cdots.   \label{Rabi-continued-fraction1}
\eeq
%Note that the ratio $T_0=\frac{{\cal R}_1^{min}}{{\cal R}_0^{min}}$ involves ${\cal R}_0^{min}$, 
%although the above continued fraction expression is obtained from the 2nd equation 
%of (\ref{Rabi-3term}), i.e the recurrence (\ref{Rabi-3term}) for $n\geq 1$. 
However, for the single-ended sequences appearing in the infinite series expansions
(\ref{Rabi-series-solution}),  the first equation (i.e. the $n=0$ part) of 
the recurrence (\ref{Rabi-3term}) requires that
\beq
 T_0=-X_0=-\frac{1}{2g}\left(E+\delta-\frac{3g^2}{\omega}+\frac{\Delta^2}
    {E-\delta+\frac{g^2}{\omega}}\right).   \label{Rabi-continued-fraction2}
\eeq
In general, (\ref{Rabi-continued-fraction1})  and (\ref{Rabi-continued-fraction2}) can
not be both satisfied for arbitrary values of the recurrence coefficients $X_n$ and $Y_n$.  
%As a result, general solutions to the recurrence (\ref{Rabi-3term}) are dominant and are usually
%generated by simple forward recursion from a given value of ${\cal R}_0$. 
Physical meaningful solutions are those that are entire in the Bargmann-Hilbert space
\cite{Schweber67}. They can be obtained if $E$
can be adjusted so that equations (\ref{Rabi-continued-fraction1}) and (\ref{Rabi-continued-fraction2})
are both satisfied. Then the resulting solution sequence ${\cal R}_n(E)$ will be purely minimal and
the power series expansions (\ref{Rabi-series-solution}) will converge in the whole complex plane.
Equating the right hand sides of (\ref{Rabi-continued-fraction1}) and (\ref{Rabi-continued-fraction2})
yields a transcendental equation for the spectrum $E$. That is,
regular energies of the driven Rabi model are given by the zeros of the transcendental function
$Q(E)=T_0+X_0$, where $T_0$ is the continued fraction in (\ref{Rabi-continued-fraction1}) and $X_0$ 
is the right hand side of (\ref{Rabi-continued-fraction2}). 
The transcendental eigenvalue equation $Q(E)=0$  may be solved for $E$ by standard root-search techniques. 






\bebb{99}

\bbit{Vedral06}
V. Vedral, Modern foundations of quantum optics, Imperial College Press, London, 2006.

\bbit{Khitrova06}
G. Khitrova,  H.M. Gibbs, M. Kira, S.W. Koch and A. Scherer, 
Nature Phys. {\bf 2}, 81 (2006).

\bbit{Thanopulos04}
I. Thanopulos, E. Paspalakis and Z. Kis, Chem Phys. Lett. {\bf 390}, 228 (2004).

\bbit{Englund07}
D. Englund et al, Nature {\bf 450}, 857 (2007).

\bbit{Niemczyk10}
T. Niemczyk et al, Nature Phys. {\bf 6}, 772 (2010).

%\bbit{Leibfried03}
%D. Leibfried, R. Blatt, C. Monroe and D. Wineland, Rev. Mod. Phys. {\bf 75}, 281 (2003).

\bbit{Schweber67}
S. Schweber, Ann. Phys. {\bf 41}, 205 (1967).

\bbit{Braak11}
D. Braak, Phys. Rev. Lett. {\bf 107}, 100401 (2011).

\bbit{Solano11}
E. Solano, Physics {\bf 4}, 68 (2011).

\bbit{Moroz12}
A. Moroz, Europhys. Lett. {\bf 100}, 60010 (2012). %\\
%Comment on ``Integrability of the Rabi model", 
%arXiv:1205.3139 [quant-ph].

\bbit{Maciejewski12}
A.J. Maciejewski, M. Przybylska and T. Stachowiak, %How to calculate spectra of Rabi and related models,
arXiv:1210.1130 [math-ph];  %\\
%Comment on ``Note on the analytic solution of the Rabi model", 
arXiv:1211.4639 [quant-ph].

\bbit{Travenec12}
I. Trav\'enec, Phys. Rev. A {\bf 85}, 043805 (2012).

\bbit{Chen12}
Q.H. Chen, C. Wang, S. He, T. Liu and K.L. Wang, Phys. Rev. A {\bf 86}, 023822 (2012).

\bbit{Moroz13}
A. Moroz, %On solvability and integrability of the Rabi model,  
Ann Phys. {\bf 338}, 319 (2013).

\bbit{Zhang13a}
Y.-Z. Zhang,  J. Phys. A {\bf 46}, 455302 (2013).

\bbit{Zhang13b}
Y.-Z. Zhang, J. Math. Phys. {\bf 54}, 102104 (2013).

\bbit{Zhong13}
H. Zhong, Q. Xie, M.T. Batchelor and C. Lee, J. Phys. A {\bf 46}, 415302 (2013).

\bbit{Moroz14}
A. Moroz, Ann Phys. {\bf 340}, 252 (2014).

\bbit{Brune87}
M. Brune,  J.M. Raimond, P. Goy, L. Davidovich and S. Haroche, 
Phys. Rev. Lett. {\bf 59}, 1899 (1987).

\bbit{Valle10}
E. del Valle,   S. Zippilli, F.P. Laussy, A. Gonzalez-Tudela, G. Morigi and C. Tejedor,
Phys. Rev. B {\bf 81}, 035302 (2012).

\bbit{Ota11}
Y. Ota,  S. Iwamoto, N. Kumagai and Y. Arakawa, 
%Spontaneous two photon emission from a single quantum dot,
 arXiv:1107.0372 [quant-ph].

\bbit{Emary02}
C. Emary and R.F. Bishop, J. Phys. A {\bf 35}, 8231 (2002).

\bbit{Gautschi67}
W. Gautschi, SIAM Review {\bf 9}, 24 (1967).

\bbit{Leaver86}
E.W. Leaver, J. Math. Phys. {\bf 27}, 1238 (1986).

\bbit{Liu92}
J.W. Liu, J. Math. Phys. {\bf 33}, 4026 (1992).

\bbit{Larson13}
J. Larson, J. Phys. B {\bf 46}, 224016 (2013).


\eebb



\end{document}

%%%%%%%%%5555555%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%666666666%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%






\bbit{Reik82}
H.G. Reik, H. Nusser and L.A. Ribeiro, J. Phys. A: Math. gen. {\bf 15}, 3431 (1982).

\bbit{Kus85}
M. Kus, J. Math. Phys. {\bf 26}, 2792 (1985).

\bbit{Lee10}
Y.-H. Lee, W.-L. Yang and Y.-Z. Zhang, J. Phys. A: Math. Theor. {\bf 43}, 185204 (2010); ibid {\bf 43}, 375211 (2010).

\bbit{Lee11}
Y.-H. Lee, J.R. Links and Y.-Z. Zhang, Nonlinearity {\bf 24}, 1975 (2011).
1304.7827.extracted_bib
\begin{thebibliography}}
\newcommand {\eebb}{\end{thebibliography}