Matrix exponential for non-commutative operator entries of matrix I would like to find the matrix exponential $e^{iHt}$ of the Hermitian matrix $H$ where
$$
H=\begin{pmatrix}
\delta& \sqrt{2}a & 0\\
\sqrt{2}a^\dagger &0& \sqrt{2}a\\
0 &\sqrt{2}a^\dagger&-\delta
\end{pmatrix}
$$
$$[a,a^\dagger]=a a^\dagger-a^\dagger a=1$$
$$[a,\delta]=[a^\dagger,\delta]=0.$$
The dagger $\dagger$ denotes complex conjugate. My thought is to make use of the definition of matrix exponential $e^{X}=\sum_{k=0}^{\infty} \frac{1}{k !} X^{k}$ and calculate $H^2$ and $H^3$ to find some pattern to exploit. For example, in the special case where $\delta=0$, we can find a diagonal matrix $D$
$$
D=\begin{pmatrix}
4a^\dagger a+6 &0&0\\
0&4a^\dagger a+2&0\\
0&0&4a^\dagger a-2
\end{pmatrix}
$$
such that $H^{2n+1}=D^{n}H$ and $H^{2n+2}=D^{n}H^2$, with which we can obtain a neat result of $e^{i H t}$. However, it seems that no such diagonal matrix exists for $\delta \neq 0$ case. Are there any other ways to calculating matrix exponential like this with non-commuting entries? Any input is much appreciated!
 A: This is really just an extended comment. I'll introduce matrices
$$S_1
=\frac{1}{\sqrt{2}}\begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{pmatrix},\quad
S_2
=\frac{1}{\sqrt{2}}\begin{pmatrix} 0 & -i & 0 \\ i & 0 & -i \\ 0 & i & 0 \end{pmatrix},\quad
S_3
=\begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1 \end{pmatrix}.$$
which satisfy the familiar angular momentum algebra $[S_j,S_k]=i\epsilon_{jkl}S_l$. A quick calculation shows
$$S_1^2+S_2^2+S_3^2=2I_3,$$
corresponding to a spin of $s=1$. Continuing in this spirit, let $S_{\pm}:=S_1\pm i S_2$. This yields
$$[S_3,S_{\pm}]=\pm S_3,\qquad [S_+,S_-]=2S_3.$$
This renders the Hamiltonian as $H=a\,S_+ + a^\dagger \,S_-+\delta\,S_3$. One useful feature of this form is that $S_\pm^3=0$, which means that powers of $H$ have some possibility of simplifying.

One additional observation. The number operator $a^\dagger a I_3$ commutes with $S_3$, so we have a basis of states $|n,m_s\rangle=|n\rangle\otimes |m_s\rangle$ for nonnegative integer $n$ and $m_s=1,0,-1$. Neither of these commutes with the Hamiltonian. However, note that the action of $aS_+$ on a state $|n,m_s\rangle$ is to reduce $n$ by $1$ and raise $m_s$ by $1$ (if possible):
$$aS_+|n,m_s\rangle = a |n\rangle \otimes S_+|m_s\rangle \propto |n-1,m_s+1\rangle$$
So the action of $aS_+$  is to move one quanta from $n$ to $m_s$ and so does not change the value of the quantum number $n+m_s$. The same goes for $a^\dagger S_-$ (moves one quanta from $m_s$ to $n$) and $\delta S_3$ (leaves both numbers unchanged). This suggests introducing the combined number operator $N=a^\dagger a I +S_3$ where $I$ is the 3-by-3 identity matrix. A quick calculation verifies that this commutes with the Hamiltonian:
\begin{align}
[N,H]
&=[a^\dagger a I+S_3,a S_++a^\dagger S_-+\delta S_3]\\
&=[a^\dagger,a] aS_+ + a^\dagger [a,a^\dagger]S_- +a[S_3,S_+]+a^\dagger [S_3,S_-]\\
&=-aS_+ +a^\dagger S_-+a S_3-a^\dagger S_-\\&=0
\end{align}
Thus we can diagonalize $H$ and $N$ simultaneously.
A: A tip of the hat to @Semiclassical 's insightful suggestion. Now that he has unravelled the group structure of H through the symmetry operator N, you may classify its eigenstates of eigenvalue n and specify them further to find eigenstates of the hamiltonian.
Recall $a^\dagger |n\rangle=\sqrt{n+1}|n+1\rangle$ and $a|n\rangle= \sqrt{n}|n-1\rangle$, so impose $H V(n)=\epsilon(n) ~ V(n)$ to fix the labile coefficients b and c for the unnormalized N eigenstate,
$$
V(n)\equiv \begin{pmatrix} \sqrt{2} b(n) |n-1\rangle \\ |n\rangle \\\sqrt{2} c(n) | n+1\rangle    \end{pmatrix}.
$$
The Hamiltonian eigenvalue condition then restricts
$$
\delta b +\sqrt{n}=\epsilon b \\
2\sqrt{n} b + 2\sqrt{n+1} c = \epsilon\\
\sqrt{n+1}-\delta c =\epsilon c, \leadsto \\
b=\frac{\sqrt n}{\epsilon -\delta}, ~~ c=\frac{\sqrt {n+1}}{\epsilon +\delta} ,\\
\epsilon (\epsilon^2 -\delta^2 -4n-2)+2\delta =0.
$$
The (readily depressed!) cubic equation may be solved for $\epsilon$, but I won't go into it... For vanishing δ, you have
$$
\epsilon =\pm \sqrt{4n+2},\\
b(n) = 1/\sqrt{4+2/n}, ~~~ c(n)= 1/\sqrt{4-2/(n+1)},
$$
so that
$V(n,t) = \exp \left( \pm it \sqrt{4n+2}\right )~~ V(n,0) $, etc. Even if $D-H^2\neq 0$, $V(n)$ is in its kernel.
I have not fully explored the non-vanishing δ solutions of the (depressed!) cubic.
If you wished to find the operator exponentiation (but why?) ,
you could write its spectral resolution in terms of these states $V(n)$, after you normalize them.
