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.