Numerical Solution of 1D Time Dependent Schrödinger Equation by Split Operator Fourier Transform (SOFT) Method¶
Authors: Dou Du, Taylor James Baird and Sara Bonella
Source code: https://github.com/osscar-org/quantum-mechanics/blob/master/notebook/quantum-mechanics/soft.ipynb
This notebook presents the split operator Fourier transform (SOFT) numerical method for solving the one-dimensional time-dependent Schrödinger equation.
Goals¶
- Understand the various steps in the SOFT algorithm and how they are translated to code.
- Familiarize yourself with the key assumptions underlying the algorithm and how their validity depends on the values of the parameters used in the simulation.
- Use the SOFT method to investigate various quantum mechanical features which appear in the dynamics of a quantum system.
- Observe the manifestation of Heisenberg's uncertainty principle in the time evolution of a nuclear wavepacket.
Background theory¶
Tasks and exercises¶
- Investigate the dependence on timestep of the stability of the dynamics (try moving the slider for $dt$ and monitor the behaviour of the various control properties of the simulation).
Solution
One may observe that as the timestep employed in the simulation is increased, the conservation of total energy of the system degrades until eventually the timestep is so large that the dynamics becomes totally unstable. The reason why this integration scheme does not conserve total energy exactly may be attributed to the non-commutativity of the split-operator propagator with the Hamiltonian. It is worth noting, however, that norm conservation is maintained even as one increases the timestep. This latter fact is due to the unitarity of the propagator in the split-operator scheme. - What dictates the maximum size of the timestep that we can use for the SOFT algorithm (consider the main assumptions/approximations that are made when formulating the propagation scheme for SOFT).
Solution
Recall that the central approximation used in the derivation of the SOFT propagation scheme is in the truncation of the Trotter product formula: $e^{A+B} = \lim\limits_{P \to \infty} (e^{\frac{A}{P}} e^{\frac{B}{P}})^{P} \approx (e^{\frac{A}{N}} e^{\frac{B}{N}})^{N}$ for $N$ sufficiently large. This approximation becomes more and more accurate the larger we make the value of $N$. In our specific case we have $e^{\frac{it}{\hbar} (\hat{T} + \hat{V} )}$ and we approximate this via the product $(e^{\frac{idt}{\hbar}\hat{T} } e^{\frac{idt}{\hbar}\hat{V} })^{N_{\text{steps}}}$ where $N_{\text{steps}}\cdot dt = t$. - Why is the use of the SOFT algorithm unfeasible for larger systems (think about how much data is needed to represent the state of the system on a computer and how many operations are required to carry out the propagation)?
Solution
In order to implement the MSOFT algorithm on a computer it is necessary to discretize the nuclear wavefunction. To do so, we must introduce a grid of points that make up the space throughout which the wavefunction extends. Say that for each dimension of the grid we use $N$ grid points. For a system in $d$ dimensions, this means that we shall require $N^d$ points to represent the wavefunction of a single nucleus. Now, if we want to instead consider a system of, say $n$ nuclei, then we find that a total number of $N^{nd}$ grid points are required. In other words, the amount of data required to represent our system scales exponentially with increasing system size. Moreover, since we must loop over each of these grid points to carry out the time evolution of our system - the number of operations required also scales exponentially. This is what renders the SOFT algorithm unsuitable for large systems. - Observe the evolution of the wavefunction in both position and momentum space. Do you notice the manifestation of Heisenberg's uncertainty principle? Can you explain it?
Solution
The top figure shows the relation between position $x$ and the modulus of the wavefunction $|\psi|$. On the other hand, the bottom figure shows the relation between momentum $p$ and the modulus of the wavefunction, $|\psi|$. Using the Harmonic potential for example, one can clearly see that when the position of the particle is known with greater certainty (a sharp peak of the $|\psi(x)|$), the $|\psi(p)|$ will have a boarder shape. This remains true for all time.
Legend¶
(How to use the interactive visualization)
Interactive figures¶
The interactive figures consist of four subplots. Clockwise from the top left we have:
- The evolution of the wavefunction in the position representation, $x$. The chosen potential is also shown in this subplot.
- The evolution of the wavefunction in the momentum representation, $p$. The initial momentum $p_0$ (along with $-p_0$) are indicated with veritcal lines in this subplot.
- Plot of the kinetic, potential and total energies vs simulation timestep (zoom to better discern the evolution of each contribution).
- Evolution of the total norm of the nuclear wavepacket.
Controls¶
Before running the simulations, users need to select a potential and set various simulation parameters. We offer three types of potentials, namely the box, Morse and Harmonic potentials. One can also tune the parameters of these potentials. Moreover, the nuclear mass, $m$, and the size of the time step, $\Delta t$ may be modified via their corresponding sliders. Afterwards, one can click the "Play" button to run the simulation. One can also pause the simulation and reselect potential and reset parameters.
By default, the figures only show the absolute value of the wavefunction, $|\psi(x)|$. However, one can also select to show the real Re$(\psi(x))$ and image Im$(\psi(x))$ parts of the wavefunction by ticking the corresponding checkboxes.