  We present the analytical solution of the Tavis-Cummings (TC) model for more
than one qubit inhomogeneously coupled to a single mode radiation field beyond
the rotating-wave approximation (RWA). The significant advantage of the
displaced oscillator basis enables us to apply the same truncation techniques
adopted in the single qubit Jaynes-Cummings (JC) model to the multiple qubits
system. The derived analytical spectrum match perfectly the exact
diagonalization numerical solutions of the inhomogeneous TC model in the
parameter regime where the qubits transition frequencies are far off-resonance
with the field frequency and the interaction strengths reach the ultra-strong
coupling regime. The two-qubit TC model is quasi-exactly solvable because part
of the spectra can be determined exactly in the homogeneous coupling case with
two identical qubits or with symmetric(asymmetric) detuning. By means of the
fidelity of quantum states we identify several nontrivial level crossing points
in the same parity subspace, which implies that homogeneous coupled two-qubit
TC model with $\omega_1=\omega_2$ or $\omega_1\pm\omega_2=2\omega_c$ is
integrable. We further explore the time evolution of the qubit's population
inversion and the entanglement behavior taking two qubits as an example. The
analytical methods provide unexpectedly accurate results in describing the
dynamics of the qubit in the present experimentally accessible coupling regime,
showing that the collapse-revival phenomena emerge, survive, and are finally
destroyed when the coupling strength increases beyond the ultra-strong coupling
regime. The suggested procedure applies readily to the multiple qubits system
such as the GHZ state entanglement evolution and quantum entanglement between a
single photon and superconducting qubits of particular experiment interest.
