recently I've been studying the toric code (a squared lattice in the context of quantum computation). I want to calculate the energy of the ground state and of all the excitations, with the respective degeneracies. This translates in calculating the eigenvalues of the hamiltonian of the system. This hamiltonian can be expressed as a matrix $H$:
\begin{equation}
    H=-A-B=-\sum_{i=1}^{N^2}A_i-\sum_{i=1}^{N^2}B_i=-\sum_{i=1}^{N^2}\bigotimes_{j=1}^{2N^2}P_{i,j}-\sum_{i=1}^{N^2}\bigotimes_{j=1}^{2N^2}Q_{i,j}
\end{equation}
where $P_{i,j}\in\{I,X\}$ and $Q_{i,j}\in\{I,Z\}$, $\forall i,j$, and where $I,X,Z$ are the usual $2\times 2$ Pauli matrices.

Moreover, the products of Pauli matrices are such that $A_i$ and $B_i$ commute altogether, having $[A_i,A_j]=[B_i,B_j]=[A_i,B_j]=[A_i,H]=[B_i,H]=0$ $\forall i,j=1,...,N^2$. (To understand how the $A_i$'s and $B_i$'s are built and why they commute, check the following image:
[![enter image description here][1]][1]


  [1]: https://i.sstatic.net/o2NJ6.png

where we can assign a number from $j=1,...,2N^2$ to each "circle" (spin), and thus notice that each $A_i$ and $B_i$, $i=1,...,N^2$, acts only on four spins, thus, only four $P_{i,j}$, resp $Q_{i,j}$, will be different from the identity. For instance, for $N=2$, an example is $P_{1,1}=X\otimes X\otimes I\otimes I\otimes X \otimes X \otimes I\otimes I$, and $Q_{1,1}=Z\otimes I\otimes Z\otimes I\otimes Z\otimes I\otimes Z\otimes I$.)

We now come back to the main goal of finding the eigenvalues of $H$. Because $A,B$ and $H$ commute and are also normal, they can be simultaneously diagonalised. Meaning that the eigenvalues of $H$, $h_i$, $\forall i=1,...,2^{2N^2}$, can be expressed as
\begin{equation}
h_i=-\sum_{j=1}^{N^2}a_{i,j}-\sum_{j=1}^{N^2}b_{i,j},\quad\forall i=1,...,2^{2N^2},
\end{equation}
where $a_{i,j}$ and $b_{i,j}$ are the different eigenvalues of $A_j$ and $B_j$, respectively.

On the other hand, if we now consider $C$, $D$ random matrices of the same dimension $n$, and thanks to the properties of the Kronecker product, we know that the eigenvalues of $M=C\otimes D$, are $eig(M)=\{c_1d_1,...,c_1d_n,c_2d_1,...,c_2d_n,...,c_nd_1,...,c_nd_n\}$, where $c_i$ and $d_j$ are the eigenvalues of $C$ and $D$, respectively. We can extrapolate this result to calculate the eigenvalues of $A$ and $B$ as all the possible combinations of products of the eigenvalues of $A_i$ and $B_i$, respectively.

With the two above results, one should be able to calculate $h_i$ as a sum of the product of the eigenvalues of the different $P_{i,j}$ and $Q_{i,j}$.

However, I'm a bit concerned about the order of the eigenvalues of $A$ and $B$. As there are many ways to diagonalize every matrix (in the sense that they can be ordered in many ways), I am not sure that the sum made to obtain $h_i$ will give me an appropriate result. So, my question is, how can I make sure that the eigenvalues obtained using this method are the proper ones?

I am approaching this in a rather "mathematical" way. However, any "physical" perspective is also welcome. My final aim is, as I've said in the beginning, to calculate the energy of the ground state and of all the excitations, with the respective degeneracies, depending on $N$.