How to find eigenvalues of a Pauli matrix resolved direct product matrix? For a $2 \times 2$ hermitian matrix one can resolve the matrix
in terms of Pauli matrices like this
\begin{align}
 H &= \begin{pmatrix}
  a & b \\
  b & -a
 \end{pmatrix} \\
 &= a \sigma_z + b \sigma_x  \\
\end{align}
Here, I've assumed $ a,b \in \mathbb{R} $, Using Pauli matrix identities one
can find the eigenvalues of $H$ are $\pm \sqrt{a^2+b^2}$
I was thinking if this process can be generalized for higher dimensional matrices.
Lets say, we have a $(2n\times2n)$ matrix and this matrix can be written in terms of
Pauli matrices
\begin{align}
 H = \sigma_{x}\otimes \mathbf{A} + \sigma_{y}\otimes \mathbf{B}
\end{align}
For example we set $n=2$
\begin{align}
 \mathbf A = \begin{pmatrix}
  e_1 & 0 \\
  0 & e_2
 \end{pmatrix} \\
 \mathbf B = \begin{pmatrix}
  b_1 & b_2 \\
  b_2 & b_4
 \end{pmatrix}
\end{align}
Can I find the eigenvalues using Pauli matrix identities as the $2\times2$ matrix, here also?
In $2\times2$ matrix one would do
\begin{align}
 H^2 = \left(
\begin{array}{cc}
 a^2+b^2 & 0 \\
 0 & a^2+b^2 \\
\end{array}
\right) \\
E^2 =\left(   
\begin{array}{cc}
 a^2+b^2 & 0 \\
 0 & a^2+b^2 \\
\end{array}
\right) 
\end{align}
Thus one obtains the eigenvalue of H as $E=\pm \sqrt{a^2+b^2}$.
When I'm doing the same thing here
\begin{align}
 H^2 &= (\sigma_{x}\otimes \mathbf{A} + \sigma_{y}\otimes \mathbf{B})^2  \\
 \implies E^2 \mathbb{I}_{2n} &= \sigma_x\sigma_x \otimes \mathbf{A}^2 + \sigma_y\sigma_y\otimes \mathbf{B}^2
 -\sigma_x\sigma_y \otimes \mathbf{A}\mathbf{B} - \sigma_y\sigma_x\otimes \mathbf{B}\mathbf{A} \\
 \implies E^2 \mathbb{I}_{2n}&=  \mathbb{I}_2 \otimes \mathbf{A}^2 + \mathbb{I}_2\otimes \mathbf{B}^2
 +i\sigma_z \otimes \mathbf{A}\mathbf{B} - i\sigma_z\otimes \mathbf{B}\mathbf{A} \\
 \implies E^2 \mathbb{I}_{2n}&=  \mathbb{I}_2 \otimes (\mathbf{A}^2 + \mathbf{B}^2) + i \sigma_z \otimes [\mathbf{A},\mathbf{B}]
\end{align}
Here $[\mathbf{A},\mathbf{B}] = \mathbf{A}\mathbf{B} - \mathbf{B}\mathbf{A}$.
Now, when I take trace of the last equation $(n =2)$
\begin{align}
4 E^2 = Tr \mathbf{A}^2 + Tr \mathbf{B}^2 
\end{align}
There should be no other terms as Pauli matrices are traceless,
so
\begin{align}
 E = \pm \frac{1}{2}\sqrt{4 + b_1^2 + 2 b_2^2 + b_4^2 + e_1^2 + e_2^2}
\end{align}
Both of the solutions are 2 fold degenerate,
But the solutions are wrong, Mathematica gives some complicated and totally different output
\begin{align}
 \sqrt{\frac{1}{2} \left(\lambda-\sqrt{\lambda^2-4 \delta}\right)} \\
 \sqrt{\frac{1}{2} \left(\lambda-\sqrt{\lambda^2+4 \delta}\right)} \\
 - \sqrt{\frac{1}{2} \left(\lambda-\sqrt{\lambda^2-4 \delta}\right)} \\
 - \sqrt{\frac{1}{2} \left(\lambda-\sqrt{\lambda^2-4 \delta}\right)} \\
\end{align}
Here $\lambda =  Tr \mathbf{A}^2 + Tr \mathbf{B}^2$, and $\delta = \text{Det}\: H$.
What I'm doing wrong here and how to get the correct expressions?

The problem is your assumption that $H^2 = E^2 I_{2n}$.
In the $2 \times 2$ case, we use the equation $H^2 = E^2 I_2$ in order to exploit the fact that any trace-zero $2 \times 2$ matrix $H$ will be such that $H^2$ is a multiple of the identity.  However, this trick no longer works in the generalized setting since it is no longer true that a trace-zero Hermitian matrix $H$ is such that $H^2$ is a multiple of the identity.
