$\Lambda_{s}$ is the strength of the interaction between nearest-neighbour qubits. $\Omega_{p} = (\Delta_{p}^{2} + \epsilon^{2})^{1/2}$ is the energy difference between the two levels of the $p^{\rm th}$ qubit, where $\Delta_{p}$ is the inherent excitation gap and $\epsilon$ is the detuning of the external magnetic flux from the flux quantum $h/(2 e)$. In our numerical computations we have used the experimentally relevant [@QuantMeta] parameter values $\Omega_{\mathrm{f}}/(2 \pi)
 = 7.78 \,{\rm   GHz}$ and $\epsilon/(2 \pi) =  4.62\,
 {\rm  GHz}$. The level separations $\Delta_{p}$ of the individual qubits have been drawn from a Gaussian distribution with a mean given by ${\langle \Delta \rangle}/(2 \pi) =  5.6\,{\rm  GHz}$ and a standard deviation $0.2 \,{\langle \Delta \rangle}$.

We have considered three cases, namely, (i) $\Lambda_{s} = \chi = 0$ (ii) $\Lambda_{s}/(2 \pi) = 1\,{\rm MHz}, \, \chi=0$ (iii) $\Lambda_{s}/(2\pi) = 
\chi/(2\pi)  = 1\, {\rm  MHz}$. In each case, $\Lambda/(2\pi) $ is varied from $- 1.2\,{\rm  MHz}$ to $ 1.3\,{\rm  MHz}$ in steps of $0.025\,{\rm  MHz}$. It is easily shown that the total number operator $$\mathcal{N}_{\mathrm{tot}}
=a^{\dagger}a + \sum_{p=1}^{M} \sigma_{p}^{+} \sigma_{p}^{-},
\label{eq:ntotTC}$$ commutes with $H_{\textsc{tc}}$. For each value of $\Lambda$ we have numerically solved for the complete set $\lbrace \ket{\psi_{M,N,k}} \rbrace$ of common eigenstates of $\mathcal{N}_{\mathrm{tot}}$ 