The operator splitting algorithm separates the Hamiltonian into a kinetic energy term and a potential energy term so that in position space potential energy is simple function and in momentum space the kinetic energy is easily calculated.

T(k) is the kinetic energy operator, and V(x) is the potential energy operator. The solution to the time evolving Schrödinger equation is:

We would like to write separate this equation into two exponentials, the T time-evolution operator times the V time evolution operator: 

.

However, because T & V are quantum mechanical operators and might not commute, we must use the following commutation relation:

[A,B]= AB-BA

In order to ensure that the operators are not switched by the commutation of multiplication rule, we must separate the time evolution operator by symmetrically decomposing the product:

Or, in our notation:

The new (correct) wave function is:

 Now we have a Schrödinger equation without the offensive second partial differential equation and a solution separated into position and momentum space! Given an initial wave function, , we can update half of the potential energy part in position space, update the kinetic energy part in momentum space, and finally update the second half of the potential energy part back in position space.  The only problem left is how to transfer from momentum space to position space and back again.  We use the Fast Fourier Transform method.

A plane wave in momentum space has momentum . Since in our algorithms, momentum is equal to k. Therefore, when we initiate this plane wave in position space and transform it to momentum space, we expect to find only one value for momentum: k. (Remember: for the FFT to maintain accuracy by the Nyquist limit, k can only range from –numpts/2 to numpts/2.) In momentum space the transformed function, , is updated by multiplying it by the kinetic energy operator,  where . The trickiest part of this is making sure the k’s from both functions are equal. The FFT algorithm moves k values from 0..numpts/2 to the bins from 0...numpts/2, and the k values from –numpts/2…-1 to the bins from numpts/2+1…numpts. Therefore, when the kinetic energy part of the time evolution operator is calculated for each k at the beginning of the program, it must be rearranged to match the order of .

Figure 1: The graph on the left is a plane wave,  , where k=1. The transform of this wave using the FFT algorithm is on the right. The k=1 bin contains the only non-zero value and will be multiplied by exp(-iT(1)dt). When the inverse transform is performed on , the plane wave will still have a single valued k, but the phase of the real and imaginary parts will move forward.  Note that the k values on the momentum space graph appear out of order. This must be accounted for when calculating the array, T(k).

                  Aside from the ordering complications of the operator splitting algorithm, the use of the FFT changes the physical meaning of the grid space where the wave function lives. Where as the Visscher and the Crank-Nicholson grid space are flat with an infinite potential at each end, the operator splitting grid space is round so that the  highest x value cell is the same as the lowest grid cell. Because of this, the odd eigenstates of the infinite square well are not able to be simulated by using the ends as infinite potentials. The looped grid space can only maintain integer multiples of full wavelengths and still remain a steady state system.

                  The time step for the program can theoretically be set to any length. For any constant potential with respect to x, the wave will remain physically meaningful. However, if the wave propagates toward some change in potential and the time step is too large, the algorithm does not account for the effects of the change in potential at intermediate time steps. For instance, if the wave is moving toward a barrier and the time step is too large, the kinetic energy part of the time evolution operator might hurtle the wave right past the barrier or directly into the barrier without calculating any change in momentum due to the increased potential energy. Therefore, a time step should be on the order of dx or smaller so that every step in the potential is accounted for.