noninteracting Hamiltonian $-\Delta J_{z}+\omega a^{\dagger }a$. The whole Hilbert space can then be decomposed into the subspaces formed by these states which can be diagonalized analytically. It is easy to write the following tri-diagonal matrix form:

$$H_{\texttt{RWA}}=\left(
\begin{array}{cccc}
\omega (n+2)+\frac{3\Delta }{2} & T_{n+1,n+2} & 0 & 0 \\
T_{n+1,n+2} & \omega (n+1)+\frac{\Delta }{2} & T_{n,n+1} & 0 \\
0 & T_{n,n+1} & \omega n-\frac{\Delta }{2} & T_{n-1,n} \\
0 & 0 & T_{n-1,n} & \omega (n-1)-\frac{3\Delta }{2}%
\end{array}%
\right) .  \label{RWA}$$

where $$\begin{aligned}
T_{n+1,n+2}&&=g\sqrt{3(n+2)}/4,T_{n,n+1}=g\sqrt{n+1}/4, \\
T_{n-1,n}&&=g\sqrt{3n}/4.\end{aligned}$$

If CRW terms $a^{\dagger }J_{+}$ and $aJ_{-}$ are included, the Hilbert space cannot be decomposed into the finite dimensional spaces, because the total excitation number $N=a^{\dagger }a+J_{z}+3/2$ is non-conserved and the subspace for different index $n$ defined above is highly correlated. So analytical solutions in this case should be highly non-trivial.

The Hamiltonian (\[ham1\]) including the CRW terms with a rotation around the $%
y$ axis by an angle $\pi/2$ can be rewritten as $$H=\Delta J_{x}+\omega a^{\dagger }a+g(a^{\dagger }+a)J_{z}.  \label{Ham}$$Introducing a unitary transformation $U=\exp \left[ \frac{g}{\omega }%
J_{z}\left( a^{\dagger }-a\right) \right] $, one can obtain the transformed Hamiltonian $H'_{SB}=H_{0}+H_{1}$, consisting of $$\begin{aligned}
  \label{H_11}
H_{0} &=&\omega a^{\dagger }a-\frac{g^{2}}{\omega }J_{z}^{2},  \label{H_0} \\
H_{1} &=&\Delta \left\{ J_{x}\cosh 