I know that there exist probability conserving explicit solvers for time-dependent Schrödinger's equation, for example, Visscher's one.

But when I tried to take into account spin and magnetic field (so that Hamiltonian matrix is no longer real-valued) in this scheme, it appeared, although remaining stable, to not conserve total probability. I had to make time steps several orders of magnitude smaller to have total probability remain somewhat constant, and even then it visibly changed after 100 time steps.

Here's what the basic algorithm looks like. Schrödinger equation $i\frac d{dt}\Psi=H\Psi$ is rewritten with $\Psi=R+iI$ as

$$\left\{\begin{align} \frac {dR}{dt}&=HI,\\ \frac {dI}{dt}&=-HR. \end{align}\right.\tag1$$

These equations are then discretized as

$$\left\{\begin{align} R\left(t+\frac12\Delta t\right)&=R\left(t-\frac12\Delta t\right)+\Delta tHI(t),\\ I\left(t+\frac12\Delta t\right)&=I\left(t-\frac12\Delta t\right)-\Delta tHR(t). \end{align}\right.\tag2$$

Probability density is defined as

$$P(x,t)=R(x,t)^2+I\left(x,t+\frac12\Delta t\right)I\left(x,t-\frac12\Delta t\right)\tag3$$

at integer $t/\Delta t$ and

$$P(x,t)=R\left(x,t+\frac12\Delta t\right)R\left(x,t-\frac12\Delta t\right)+I(x,t)^2\tag4$$

at half-integer $t/\Delta t$ and conserved provided that $H$ matrix (where Laplacian is supposed to be a finite-difference or something similar) is real-valued.

My change was to make $H$ complex to allow for Pauli matrix $\sigma_y$ and magnetic field term $\propto i\partial_y A(\vec r)$. Now, with $H=H_r+iH_i,$ equations $(1)$ look like

$$\left\{\begin{align} \frac {dR}{dt}&=H_rI+H_iR,\\ \frac {dI}{dt}&=H_iI-H_rR. \end{align}\right.\tag5$$

After this change probability density defined by $(3)$ and $(4)$ is no longer conserved. Apparently, Visscher's algorithm is not extensible to complex Hamiltonians, so I need to find another one.

So, are there any explicit probability conserving solvers suitable for Pauli equation?

**NOTE**

There appears to be many papers related to Dirac equation, but I would like to solve the Pauli equation *itself*, not just Pauli *regime* of Dirac equation (e.g. input large potentials which would lead to Klein paradox in Dirac equation but wouldn't in Pauli equation).