  The Dirac equation is an important model in relativistic quantum mechanics.
In the semi-classical regime $\epsilon\ll1$, even a spatially spectrally
accurate time splitting method \cite{HuJi:05} requires the mesh size to be
$O(\epsilon)$, which makes the direct simulation extremely expensive. In this
paper, we present the Gaussian beam method for the Dirac equation. With the
help of an eigenvalue decomposition, the Gaussian beams can be independently
evolved along each eigenspace and summed to construct an approximate solution
of the Dirac equation. Moreover, the proposed Eulerian Gaussian beam keeps the
advantages of constructing the Hessian matrices by simply using level set
functions' derivatives. Finally, several numerical examples show the efficiency
and accuracy of the method.
