Simulation of Time-Dependent Schrödinger Equation in the Position and Momentum Domains

The paper outlines the development of a new, spectral method of simulating the Schrödinger equation in the momentum domain. The kinetic energy operator is approximated in the momentum domain by exploiting the derivative property of the Fourier transform. These results are com-pared to a position-domain simulation generated by a fourth-order, centered-space, finite-differ-ence formula. The time derivative is approximated by a four-step predictor-corrector in both domains. Free-particle and square-well simulations of the time-dependent Schrödinger equation are run in both domains to demonstrate agreement between the new, spectral methods and established, finite-difference methods. The spectral methods are shown to be accurate and precise.


Introduction
This paper outlines the development of simulations of the time-dependent Schrödinger equation produced in both position and momentum domains. In the position domain this is the x-y plane. In the momentum domain this is the k x -k y plane, as it is the Fourier transform of the position domain [1]. The simulations demonstrate the accuracy of the spectral methods used in the momentum domain.
All simulations are advanced in time using a four-step predictor-corrector method. The predictor-corrector can be applied independently in both position and momentum domains to step the simulation forward in time. The predictor-corrector is generated using Lagrange polynomials, outlined by [2] and [3]. The predictor formula found here is shown to be consistent with established Adams-Bashforth formulas [3].
The position-domain approximation of the kinetic energy operator is derived using Lagrange polynomials and consistent with results from [4]. In the position domain the approximation to the kinetic energy operator is fourth-order accurate. In the momentum domain, the kinetic energy operator approximation is global-order accurate because it relies on the derivative property of the Fourier transform [5]. The software written to generate these simulations uses the Fastest Fourier Transform in the West (FFTW) to transform between position and momentum domains [6]. Simulating the time-dependent Schrödinger equation in the momentum domain achieves higher orders of spatial accuracy. The performance and precision of momentum-domain simulations is comparable to position-domain simulations.
Given an initial state 0 y at 0 t = , the four-step predictor-corrector requires the creation of n j n y y j t − = − ∆ for 1, , 4 j =  in order to compute the first predictor-corrector time-step. A simple backwards Euler method, outlined in [4], [7]- [12], is used to generate the wave function at these early time-steps. Each of these early states for 1, , 4 j =  is re-normalized after their creation to ensure minimum initial error. The first simulation is a free particle with no imposed boundary conditions, when the Hamiltonian consists only of the kinetic energy operator. This simulation demonstrates the difference in boundary conditions of each domain. In the position domain, this is equivalent to an infinite square-well potential, or particle-in-a-box. When the wave function reaches one boundary, it is reflected back. In the momentum domain, this is equivalent to periodic boundary conditions. When the wave function disappears into one boundary it will reappear in the opposite boundary, travelling in the same direction. This is to establish a relative performance benchmark when only the kinetic energy operator is applied.
Second, a finite square-well potential of 100 eV is imposed in both domains. This simulation demonstrates the computational burden associated with imposing the same initial and boundary conditions in both domains. Application of the potential function in the position domain is carried out by entry-wise multiplication of the wave function and potential function lattices. In the momentum domain, this operation is equivalent to convolution. Rather than carry out this time-consuming operation, the wave function in the momentum domain is transformed back to the position domain at every time-step in the simulation in order to apply the potential function. The kinetic energy operator is applied when the wave function has been transformed forward into the momentum domain.
Each of the simulations begins with an initial state of a two-dimensional wave packet with a Gaussian envelope. The simulations are stepped forward in time and the complex-valued wave function components and densities, as well as some expectation values, are captured incrementally. The wave function components and densities are converted to image format and animated [13].

Methods
The following subsections outline the numerical methods used to generate solutions to the Schrödinger equation. (1)

Time Derivative
The same four-step predictor-corrector method is used to step the position and momentum domain simulations forward in time. The predictor and corrector start with the basic form of the differential equation

Predictor
The predictor uses the integral form of the differential equation.
The function f is approximated by Lagrange polynomials [2]. Once the polynomial approximation to f has been substituted, the integral is straightforward to compute. This yields Adams-Bashforth coefficients which calculate the predicted value n y  [3].

Corrector
The corrector uses the original form of the differential equation.

( )
n n y f y y f t For the corrector, y is approximated by Lagrange polynomials [2]. Once the polynomial approximation to y has been substituted, the first derivative is straightforward to compute. This yields the following coefficients which calculate the value n f .
The predicted value n y  is substituted for n y in the function f to get n f  .

Application to the Schrödinger Equation
In the position domain, the Schrödinger equation is written as follows using operator notation as shorthand. The Hamiltonian operator Ĥ is written with a P superscript to denote the position domain and with superscript M to denote the momentum domain.
At every time step, the predictor calculation is carried out.
That result is plugged in to the corrector calculation to advance the simulation forward one time-step.
The following sections outline the development of the Hamiltonian operator ˆP H in the position domain and ˆM H in the momentum domain.

Free Particle, Kinetic Energy Operator
The first simulations in the position and momentum domains assume free particle conditions. The particle is given the mass of an electron. Only the kinetic energy operator applies in the Hamiltonian.
Before the simulations are started, the position-and momentum-domain lattices must be constructed. This requires fixing the position-domain step sizes x ∆ and y ∆ as well as the number of columns x N and number of rows y N . It is also helpful to define an origin point, 2 The reversal of the y-direction accounts for the fact that computer storage increments the row index as the row moves down.
The momentum domain lattices are constructed according to the relationship 1 . The high frequencies have been shifted into the negative frequencies. Use of the FFTW library requires applying a phase-shift to the position domain before transforming into the momentum domain if negative frequencies are used instead of high frequencies [6].

Position Domain
The approximation to the kinetic energy operator in the position domain was generated using Lagrange polynomials. This was accomplished by approximating the second derivative in one dimension, as the same formula can be applied to all dimensions. This is a centered-space formula accurate to fourth order, requiring five sample points. In terms of the generalized coordinate q , the sample points and a fixed q ∆ describe the set of sample points centered around 2 q . For a function ( ) , the polynomial approximation to f is substituted and the second derivative is calculated. This procedure yields the following approximation to the Laplacian operator: 16 30 The real and imaginary parts of the wave function R iI Ψ = + are calculated independently, yielding a pair of coupled equations. 2 Denote spatial sample points with subscripts: ( )    I  I  I  I  I  m x

Momentum Domain
The approximation to the kinetic energy operator in the momentum domain was generated using the transform of the derivative operator. Because the momentum domain is the Fourier transform of the position domain, the derivative operator is transformed as well. In terms of the generalized position-domain coordinate q , let the function ( ) F s be the Fourier transform of the function ( ) and s is the generalized momentum-domain coordinate.
The initial position-domain state Ψ is transformed forward into the momentum-domain state Φ , where ( )( )

Square Well, Potential Energy Operator
A square-well potential was chosen to test the effectiveness and relative performance of the simulations of both domains, when the particle interacts with an electrostatic potential. The particle is again given the mass of an electron. For the purposes of this demonstration, the chosen potential must be high enough to reflect most of the particle off the potential step, back into the region where the potential is zero.
The lattice describing the potential must have the number of columns x N and number of rows y N . Before the simulations are started, it is helpful to choose a well boundary index constant step j and potential step size

Position Domain
The application of the potential operator is straightforward in the position domain. The lattices representing the real and imaginary parts of the wave-function Ψ are multiplied entry-wise by the lattice representing the potential function. The simulation can be stepped forward in time once the approximation to the Hamiltonian has been substituted into the predictor corrector formula.

Momentum Domain
In the momentum domain, application of the potential operator transforms to the convolution operation,  and * represents convolution [5]. Rather than carry out this computationally expensive operation, the operator is applied in the position domain. This requires transforming back and forth between the position and momentum domains at every time-step. Use of the predictor-corrector complicates this procedure somewhat, because the potential operator must be applied at the predictor step and the corrector step.
For simplicity and readability the procedure below does not write the state functions Ψ and Φ decomposed into their real and imaginary components.
First, the potential operator is applied in the position domain and transformed to the momentum domain.
The kinetic energy operator is applied to Φ and added to * Φ  to produce the predictor Hamiltonian ˆM pred H .

( )( )
Following the example of Equation (9) in Section 2.1.3, the corrector formula is applied to produce the corrected value Φ .

Results
All simulations are run under the same constraints and initial conditions to illustrate similarities and differences in the particle's position over time. The free particle simulations show differences in position due to the differences in boundary conditions between the two domains, despite having the same initial conditions. The square well simulations will show strong agreement in position due to the boundary conditions and initial conditions being the same. Regarding precision, all simulations were shown to be accurate to 8 decimal places. This value is stable for the duration of the simulations and was measured by finding the difference between the current state's density function and 1.
The physical constants of electron mass m , electron charge q and Planck's constant  are given in standard units. This helps scale the simulations to real-world dimensions, although the simulated particles are much larger than real electrons to show detail. For the square-well simulations, the electrostatic potential is set to +100.0 eV, although this is also expressed in standard units. The well boundary is established at 20 index units inside the lattice boundaries.
The row and column sizes are set to 256 × seconds. Each simulation is run for 100,000 time steps. The lifetime of each simulation is 0.1 picoseconds. The components, densities and expectation values are measured and recorded every 500 time steps. This increment is referred to as a "frame" in the graphs below and each frame is equal to 16 5.0 10 − × seconds. All initial condition parameters are given in index units and the actual parameters are multiplied by the appropriate lattice step size. The initial conditions set the particle on the +x-axis at index +70.0, relative to the origin point and close to the vertical boundary at the end of the +x-axis. The particle is given positive momentum, which will propel the particle into the boundary. The particle's spatial wavelength λ is set to 10.0 index units. The particle's predicted velocity can be found by the conservation of momentum.
The initial velocity of the particle in all simulations is 5 3.64 10 × meters per second in the +x-direction. The initial momentum x k is therefore 9 π 10 × kilogram-meters per second. In units normalized by 2π , as shown in the momentum domain simulations, the initial momentum is

Free Particle
Based on the initial conditions and predicted velocity, the particle is expected to reach the boundary at time 14 3.19 10 t − = × seconds, corresponding to the 64 th frame of the simulation. The position-domain particle reflects off the boundary while the momentum-domain particle travels through the boundary and reappears in the opposite boundary, travelling in the same direction at the same velocity. Figure 1 charts the relative position along the x-axis of the position-and momentum-domain particles. The position is given by the expectation value X . A marker has been placed at the 64 th frame to show where the paths are expected to diverge. Four free-particle animations are produced for the position domain. One overhead view and one crosssectional view are produced for the finite-difference free-particle simulation and one overhead view and one cross-sectional view are produced for the spectral free-particle simulation. Figure 2 shows the cross-sectional view of the free particle's components and density produced by finitedifference methods. The cross-section is along the x-axis in the position domain. Figure 3 shows the overhead view of the free particle's position-domain density produced by finite-difference methods. The particle diffracts as it reflects off the boundary. Figure 4 shows the cross-sectional view of the free particle's components and density produced by spectral methods. The cross-section is along the x-axis in the position domain.    . Cross-sectional animation still at 64 th frame of free particle simulation produced by spectral methods. Figure 5 shows the overhead view of the free particle's position-domain density produced by spectral methods. The particle continues to diffuse in space as it pass through the boundary.

Square Well
Based on the initial conditions and predicted velocity, the particle is expected to reach the boundary at time   . Comparing X of bound particle in square well produced by finitedifference and spectral methods.

Position Domain
Four square-well animations are produced for the position domain. One overhead view and one cross-sectional view are produced for the finite-difference square-well simulation and one overhead view and one crosssectional view are produced for the spectral square-well simulation. Figure 7 shows the cross-sectional view of the bound particle's components and density produced by finitedifference methods. The cross-section is along the x-axis in the position domain. Figure 8 shows the overhead view of the bound particle's position-domain density produced by finitedifference methods. The particle diffracts as it reflects off the potential wall. Figure 9 shows the cross-sectional view of the bound particle's components and density produced by spectral methods. The cross-section is along the x-axis in the position domain. Figure 10 shows the overhead view of the bound particle's position-domain density produced by spectral methods. The particle diffracts as it reflects off the potential wall.

Momentum Domain
One additional animation is produced for the momentum-domain, square-well simulation that shows the crosssectional view of the momentum density function. The cross-section is taken along the x k axis. Figure 11 shows a cross-sectional view of the density function in the momentum domain at the 42 nd frame. As the particle interacts with the electrostatic potential, the particle is reflected and reverses direction. In the momentum domain, this is indicated by the density function disappearing from the positive axis and reappearing on the negative axis. Figure 12 illustrates this reversal of direction over time by measuring the expectation value x K . A marker has been placed at the 42 nd frame indicating when the particle interacts with the boundary and reverses direction.

Discussion
Momentum-domain simulations of the time-dependent Schrödinger equation provide precise and accurate results; however, the application of these techniques is not limited to the Schrödinger equation. The methods described in this paper are also suitable to simulate the heat equation , where f describes temperature in space and time and α is the thermal diffusivity. The spectral methods described here may be applied to any parabolic differential equation. The spectral methods also transform the Hartree-Fock operation in a many-body problem from convolution in the position domain to entry-wise multiplication in the momentum domain, although this application is not explored here due to resource constraints.
These simulations were produced on single-core, AMD Athlon X2 processor. The spectral methods demonstrated faster performance for the free-particle simulation, while the finite difference methods demonstrated faster performance for the square-well simulation. None of the simulations employed parallel computing techniques due to the limitations of the hardware. Multiple cores would allow multiple Fourier transforms to be calculated at the same time. Because the spectral-method, square-well simulation requires multiple Fourier transforms at every time step, introducing parallel computing techniques would increase performance substantially.