Let's imagine that our Hamiltonian acts as follows
\(\mathcal{H}: \mathcal{H}_1 \otimes \mathcal{H}_2 \to \mathcal{H}_1 \otimes \mathcal{H}_2\)
For simplicity, let's suppose that
\(\mathcal{H} = (a^\dag b + b^\dag a), \)
where \(a, \,\, a^\dag\) \(\)are annihilation and creation operators acting in \(\mathcal{H_1}\) (similarly for b in \(\mathcal{H_2}\)).
We can calculate matrix element
\(<l,p|\mathcal{H}|n, k> = \sqrt{n+1} \sqrt{k} \delta_{l, n+1} \delta_{p,k-1} + \sqrt{n} \sqrt{k+1} \delta_{l,n-1} \delta_{p,k+1}\).
It (matrix of the operator) gives us, formally, 4-dimensional array. Is there a convenient way, how to rewrite it in form of square matrix? Because I'm interested in eigenvectors and in this form it's very uncomfortable to work with the 4-dimensional array.