1%
\end{array}%
\right) ,  \label{eigenstates}\end{aligned}$$where $$K_{i,n}=\{%
\begin{array}{c}
\frac{1}{\sqrt{3}B_{n}}[-\frac{2g^{2}}{\omega }-(-1)^{i}B_{n}+4\chi
_{i,n}],(i=1,2) \\
\frac{1}{\sqrt{3}B_{n}}[-\frac{2g^{2}}{\omega }-(-1)^{i}B_{n}-4\chi
_{i-2,n}],(i=3,4)%
\end{array}%
,  \label{ki}$$and $$\chi _{i,n}=\sqrt{\frac{g^{4}}{4\omega ^{2}}+(-1)^{i}\frac{g^{2}}{4\omega }%
B_{n}+\frac{B_{n}^{2}}{4}}(i=1,2)$$ with $B_{n}=\Delta G_{0}(n)$. Interestingly, the zeroth-order approximation is similar to the adiabatic approximation in the two-qubit system [@agarwal], where the transition between different manifolds is not considered, and the $n$th state is only limited to the same $n$-th manifold.

The validity of the zeroth-order approximation is restricted to the large detuning regime $\Delta/\omega\ll1$. In the zero detuning limiting, $\Delta =0$, within the same manifold $n$, $|\pm \frac{3}{%
2}\rangle |n\rangle $ and $|\pm \frac{1}{2}\rangle |n\rangle $ are nearly degenerate. For a large detuning $\Delta/\omega\ll1 $, it is reasonable to consider the qubit states with the same $n$ manifold coupled by the interaction. Especially, for a strong coupling strength $g/\omega\gg\Delta/\omega$, the diagonal terms of the approximated Hamiltonian in Eq.( \[h0\]) play a more dominant role than the off-diagonal terms dependent on $\Delta$. And the high order terms in Eq.( \[H\_11\]) still can be neglected even in the strong coupling regimes. Hence, the zeroth-order approximation is expected to work well from weak to strong coupling regimes for the large detuning case $\Delta/\omega\ll1$.

The zeroth-order energy spectrum is plotted in Fig. [energy 