How to decouple this system of 2nd order partial differential equations? Solving a problem I found this system of equations:
$$  (\partial _t ^2 + \partial _x ^2 + \partial_y ^2 + \partial_z ^2)a_1(\vec{x},t) - 4gB(x\partial_y - y\partial_x)a_{\color{Red}2}(\vec{x},t) - 4g^2B^2(x^2 + y^2)a_1 - \lambda a_1 = 0 $$
$$ (\partial _t ^2 + \partial _x ^2 + \partial_y ^2 + \partial_z ^2)a_2(\vec{x},t) + 4gB(x\partial_y - y\partial_x)a_{\color{Red}1}(\vec{x},t) - 4g^2B^2(x^2 + y^2)a_2 - \lambda a_2 = 0  $$
taking into account that  $g,B,\lambda=$ constants,
Any ideas to separate these equations in terms of $a_1$ and $a_2$?
Thanks a lot!

Here's a linear algebra approach. I'll rearrange this slightly to make the next steps a bit clearer. I'm going to use $\Delta$ to represent $\partial_t^2 + \partial_x^2 + \partial_y^2 + \partial_z^2$. If it were $-\Delta_t^2$, I'd use the standard $\square$ d'Alembert notation instead. I'll also write $x\partial_y - y\partial_x$ as $iL_z$, taking cues from quantum mechanics, where $L$ is the angular momentum operator. (I suspect this is a relativistic quantum particle in an external magnetic field anyway..)
\begin{align}
(\Delta - 4g^2 B^2(x^2+y^2) - \lambda) a_1 - 4igB L_z a_2 &= 0 \\
(\Delta - 4g^2 B^2(x^2+y^2) - \lambda) a_2 + 4igB L_z a_1 &= 0
\end{align}
Writing this as a matrix equation, this would be
$$ \begin{pmatrix} \Delta - 4g^2 B^2 (x^2+y^2) & -4igB L_z \\ 4igB L_z & \Delta - 4g^2 B^2(x^2+y^2) \end{pmatrix} \begin{pmatrix} a_1 \\ a_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} $$
This matrix is of the form
$$ \begin{pmatrix} a & -b \\ b & a \end{pmatrix}. $$
A quick check shows that this matrix is diagonalized by
$$ P = \frac{1}{\sqrt{2}} \begin{pmatrix} i & 1 \\ -i & 1 \end{pmatrix}, \qquad P^{-1} = \frac{1}{\sqrt{2}} \begin{pmatrix} -i & i \\ 1 & 1 \end{pmatrix} $$
so we can turn the above matrix equation into..
$$ \frac{1}{\sqrt{2}} \begin{pmatrix} i & 1 \\ -i & 1 \end{pmatrix} \begin{pmatrix} \Delta - 4g^2 B^2 (x^2+y^2) & -4igB L_z \\ 4igB L_z & \Delta - 4g^2 B^2(x^2+y^2) \end{pmatrix} \begin{pmatrix} a_1 \\ a_2 \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} i & 1 \\ -i & 1 \end{pmatrix} \begin{pmatrix} 0 \\ 0 \end{pmatrix} $$
simply by multiplying on the left by $P$. Multiplying by $P^{-1}P$ on the left of $\begin{pmatrix} a_1 \\ a_2 \end{pmatrix}$, we get
$$ \begin{pmatrix} \Delta - 4g^2 B^2 (x^2+y^2) + 4gB L_z & 0 \\ 0 & \Delta - 4g^2 B^2 (x^2+y^2) - 4gB L_z\end{pmatrix} \begin{pmatrix} ia_1 + a_2 \\ -ia_1 + a_2\end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} $$
This matches with what you have. The strength of this approach is it gives you a prescription to go by. Decoupling a system of linear ODEs/PDEs is tantamount to diagonalizing the matrix-of-operators system.
