I am trying to evaluate $$\exp(-ic \alpha pt /\hslash - i\beta mc^2t/\hslash),$$ where $$ \alpha = \begin{pmatrix} 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \end{pmatrix}, \quad \beta = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 \end{pmatrix}. $$ The most basic property of these matrices is, besides $\beta$ being diagonal, is $\alpha^2 = \beta^2 = \mathbb{I}_4$. We also have $\alpha \beta + \beta \alpha = 0$.
This matrix exponential arises in the study of time propagation by the 1-D free Dirac equation in relativistic quantum mechanics, which reads $$ i\hslash \partial_t \psi(p,t) = (c\alpha p + \beta mc^2) \psi(p,0),$$ where $\psi$ denotes the spinor in momentum space. This leads to the matrix exponential given above such that $$\psi(p,t) = \exp(-ic \alpha pt /\hslash - i\beta mc^2t/\hslash) \psi(p,0).$$
In the appendix of a publication [1], the solution is directly given by (with $\hslash = 1$) $$\left( \mathbb{I}_4 cos(Et) - i\frac{c\alpha p + \beta mc^2}{E} sin(Et)\right), \quad E = \sqrt{p^2c^2 + m^2 c^4},$$ which looks fairly simple, and is reminiscent of Euler's formula. However, I was still not able to derive the solution with the few methods I know. For instance, the case $X = A + N$ ,where $A$ is diagonal, $N$ is nilpotent and $AN = NA$ does not apply here because of $\alpha$, which is not nilpotent, and $\alpha$ and $\beta$ do not commute, either.
Can you see it at a glance?
[1] : see equations B.3 and C.3 on pages 39 and 41 in https://arxiv.org/abs/1107.4650 .
Taking $a=-cpt/\hbar$ and $b=-mc^2t/\hbar$ for convenience, note that $\exp(ia\alpha+ib\beta)=\sum\limits_{n=0}^{\infty}\frac{1}{n!}(ia\alpha+ib\beta)^n$. But from the fact that $\alpha^2=\beta^2$ is the identity, and $\alpha\beta+\beta\alpha=0$, we have that $(ia\alpha+ib\beta)^{2k}=(-1)^k(a^2+b^2)^k I$ (where $I$ is the identity matrix), and $(ia\alpha+ib\beta)^{2k+1}=(-1)^k i(a^2+b^2)^k(a\alpha+b\beta)$. So, gathering together the even terms yields: $\cosh(\sqrt{-a^2-b^2})I$ and gathering together the odd terms yields: $\frac{i\sinh(\sqrt{-a^2-b^2})}{\sqrt{-a^2-b^2}}(a\alpha+b\beta)$. Note that $a$ and $b$ are real, so we can pull a factor of $i$ from each square root. Recalling that $\cosh(ix)=\cos(x)$ and $\sinh(ix)=i\sin(x)$ and adding the two pieces together then yields the result.