A Time Splitting Method for the Three-Dimensional Linear Pauli Equation

The Potential Step

Step (i) of Algorithm 2.1 consists in finding, for all grid points x 𝒋 , the solution to the initial value problem

where w 𝒋 = ( w 𝒋 ⁢ 1 , w 𝒋 ⁢ 2 ) and

Then the solution of the potential step is given by U n ⁣ * ⁢ ( x 𝒋 ) = e Δ ⁢ t ⁢ 𝔅 ⁡ U n ⁢ ( x 𝒋 ) := w 𝒋 ⁢ ( Δ ⁢ t ) . For time-independent magnetic field and potentials, an analytical solution is available for all time-steps outside of the solution loop, whereas for time-dependent data the solution has to be re-computed in each time-step. In the latter case the solution can be obtained with any highly efficient ODE solver.

The Kinetic Step

In Step (ii) of Algorithm 2.1, one has to solve the initial boundary value problem

which consists of nothing but two decoupled free Schrödinger equations with periodic boundary conditions for w = ( w 1 , w 2 ) . Then the solution of the kinetic step is given by U n ⁣ * * = e Δ ⁢ t ⁢ 𝔄 ⁡ U n ⁣ * := w ⁢ ( t n + 1 ) . Hence, we can use any of the available highly efficient methods for the free Schrödinger equation. In light of the advection step, a good way is to solve the equation in Fourier space using FFT. In particular, as U n ⁣ * * = e Δ ⁢ t ⁢ 𝔄 ⁡ U n ⁣ * , we find that

Then, instead of performing an iFFT to move back to physical space, we can directly pass the Fourier space data to the next step.

The Advection Step

This substep is the most subtle step of the operator splitting method, as standard methods are usually stable only under restrictive CFL-type conditions that prevent the use of large time-step sizes. However, since it is analogous to the magnetic Schrödinger equation case, we can adapt methods in [9, 5, 12] for the 2-spinor case. We opt for the method of characteristics to solve this equation combined with Fourier interpolation. Step (iii) of Algorithm 2.1 consists of the solution of

For each of the two components of w = ( w 1 , w 2 ) and each 𝒋 , the characteristic z 𝒋 ⁢ ( ⋅ ) through x 𝒋 solves the problem

with end value prescribed at t = t n + 1 . Solving the above characteristic equation for each grid point x 𝒋 would yield the sought approximation U n ⁣ * ⁣ * * ⁢ ( x 𝒋 ) via

However, the point z 𝒋 ⁢ ( t n ) is not a grid point in general, so we do not have immediate access to the value U n ⁣ * * ⁢ ( z 𝒋 ⁢ ( t n ) ) . We need to use an interpolation method to approximate U n ⁣ * * ⁢ ( z 𝒋 ⁢ ( t n ) ) based on the knowledge of U n ⁣ * * at grid points. Since the previous step passes Fourier data to the advection step, it is natural to use Fourier interpolation to accomplish this. Following [5, Section 5.1], we evaluate a Fourier interpolation at x = z 𝒋 ⁢ ( t n ) , where the coefficients { U ^ k 1 , k 2 , k 3 n ⁣ * * } are known from step (ii) of Algorithm 2.1. In general, further choices are required to make such a trigonometric interpolation unique in a sensible way (see, e.g., [10]) but we omit discussion of this here – minimally oscillatory interpolations are usually to be preferred. Besides this uniform trigonometric method, one could employ other methods for the interpolation, e.g., the computationally more efficient non-uniform NUFFT-based approaches as in [5, Section 5.3] and [12, Section 2.2].

The Coupling Step

The coupling step contains the off-diagonal components of the Pauli equation. Step (iv) of Algorithm 2.1 consists in finding, for all grid points x 𝒋 , the solution of the following initial value problem:

where w 𝒋 = ( w 𝒋 ⁢ 1 , w 𝒋 ⁢ 2 ) and

Then the solution of the coupling step, which is also the approximation U n + 1 ≈ u ⁢ ( t n + 1 ) , is given by

Unlike the previous steps, this is a coupled system of ODEs, which may be treated with appropriate highly efficient solvers. An analytic solution to this ODE is readily available in each time step, and as with the potential step (step (i) of Algorithm 2.1), for the case of time-independent potentials the solution operator may in fact be pre-computed for all considered time-steps outside of the solution loop.