(annihilation) operator $%
a^{\dagger }\left( a\right) $ appears in the GRWA, so the qubits states with different oscillator number $n$, $n\pm1$ and $n+2$ are coupled with each other. In the basis of $|-\frac{3}{2}\rangle |n+2\rangle $, $|-\frac{1}{2}\rangle
|n+1\rangle $, $|\frac{1}{2}\rangle |n\rangle $ and $|\frac{3}{2}\rangle
|n-1\rangle $ ($n>0$), the Hamiltonian $H_{\texttt{GRWA}}$ can be written in the matrix form as

$$H_{\texttt{GRWA}}=\left(
\begin{array}{cccc}
\omega (n+2)+\mu _{1}(n+2) & \Delta R_{n+1,n+2}^{\prime } & 0 & 0 \\
\Delta R_{n+1,n+2}^{\prime } & \omega (n+1)+\mu _{2}(n+1) & \Delta R_{n,n+1}^{\prime }
& 0 \\
0 & \Delta R_{n,n+1}^{\prime } & \omega n+\mu _{3}(n) & \Delta R_{n-1,n}^{\prime } \\
0 & 0 & \Delta R_{n-1,n}^{\prime } & \omega (n-1)+\mu _{4}(n-1)%
\end{array}%
\right) ,$$

with $R_{n+1,n+2}^{\prime }=\frac{-\sqrt{3}K_{2}+K_{1}(\sqrt{3}%
+2K_{2})}{C_{1}C_{2}}R_{n+1,n+2}\sqrt{n+2}$, $R_{n,n+1}^{\prime }=\frac{-%
\sqrt{3}K_{3}+K_{2}(\sqrt{3}-2K_{3})}{C_{2}C_{3}}R_{n,n+1}\sqrt{n+1}$ and $%
R_{n-1,n}^{\prime }=\frac{-\sqrt{3}K_{4}+K_{3}(\sqrt{3}+2K_{4})}{C_{3}C_{4}}%
R_{n-1,n}\sqrt{n}$.

To this end, the GRWA can be also performed analytically without more efforts than those in the original Hamiltonian $H_{\texttt{RWA}}$ in Eq.(\[RWA\]). The displaced oscillator states $|n\rangle _{m}$, $|n\pm
1\rangle _{m}$ and $|n+2\rangle _{m}$ depend upon the Dicke state $%
|j,m\rangle $, and are definitely different from both the RWA ones and the zeroth-order approximations where only the state $|n\rangle_{m}$ is considered. Hence, as $\Delta/\omega$ increases, the first-order correction provides an efficient, yet accurate analytical 