Do you need the matrix exponential, or the action of it on a vector? Since it is a Hamiltonian from quantum mechanics that you’re dealing with, you can just use one of the many approaches to propagate wavefunctions that are available in the literature. The Hamiltonian you have is Hermitian and purely real, the complexity comes about from the -im*t in the propagator, but this does not change any properties of the underlying Hamiltonian, except rotating the spectrum by -pi/2, i.e. all the real eigenvalues are now purely imaginary.
The simplest, but still stable (A-stable in fact) method that I know of is Crank–Nicolson:
which has a local error of t^3 and a global error of t^2. It amounts to multiplication by tridiagonal matrix, and solving by another [whose factorization you can easily precompute, factorize(I + im*H*t/2)]. This is a Padé approximant to the exponential function, as @stevengj mentioned.
There are also Krylov methods that can help you compute the action, see e.g. ExponentialUtilities.jl, but beware that they struggle when your grid spacing becomes too fine.