Quantum Harmonic Oscillator Wavepacket Simulation

This project is a numerical simulation of Gaussian wavepacket dynamics in a one-dimensional quantum harmonic oscillator.

The time-dependent Schrödinger equation is solved using the Split-Step Fourier method (SSFM). The purpose of the project is to study how a localized quantum wavepacket evolves in time inside a harmonic potential, how the expectation value of position changes, and how well the numerical method preserves probability.

---

Overview

In quantum mechanics, the harmonic oscillator is one of the most important model systems. It appears in many areas of physics and provides a useful test case for numerical methods.

In this project, a Gaussian wavepacket is placed away from equilibrium inside a one-dimensional harmonic oscillator potential. The wavefunction is then propagated in time using the Split-Step Fourier method. During the simulation, the code tracks:

• the probability density
• the expectation value of position
• the normalization of the wavefunction
• the normalization error over time

This project was made as a computational study and implementation of wavepacket dynamics in a standard quantum system.

---

Physical Model

The system is described by the time-dependent Schrödinger equation:

$$i\hbar \frac{\partial \psi(x,t)}{\partial t} = \left[ -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x) \right]\psi(x,t)$$

For the quantum harmonic oscillator, the potential is

$$V(x) = \frac{1}{2} m \omega^2 x^2$$

In this code, dimensionless units are used:

• $\hbar = 1$
• $m = 1$
• $\omega = 1$

With these choices, the potential becomes

$$V(x) = \frac{1}{2} x^2$$

and the Schrödinger equation simplifies to

$$i \frac{\partial \psi(x,t)}{\partial t} = \left[ -\frac{1}{2}\frac{\partial^2}{\partial x^2} + \frac{1}{2}x^2 \right]\psi(x,t)$$

---

Initial Wavepacket

The simulation starts with a Gaussian wavepacket of the form

$$\psi(x,0) = \exp\left( -\frac{(x-x_0)^2}{2\sigma^2} \right) \exp(i k_0 x)$$

where:

• $x_0$ is the initial position of the packet
• $\sigma$ is the width of the packet
• $k_0$ is the initial momentum

In the current simulation, the initial parameters are:

• $x_0 = -3.0$
• $\sigma = 0.7$
• $k_0 = 0.0$

This means the packet starts displaced from the center of the harmonic potential and then oscillates in time.

---

Numerical Method

The simulation uses the Split-Step Fourier method.

The Hamiltonian can be separated into two parts:

$$\hat{H} = \hat{T} + \hat{V}$$

where

$$\hat{T} = -\frac{1}{2}\frac{\partial^2}{\partial x^2}$$

is the kinetic energy operator, and

$$\hat{V} = V(x)$$

is the potential energy operator.

The main idea of the Split-Step Fourier method is:

• apply the potential evolution in position space
• transform the wavefunction into momentum space using the Fast Fourier Transform (FFT)
• apply the kinetic evolution in momentum space
• transform back to position space
• apply the potential evolution again

For one time step $\Delta t$, the evolution is approximated as

$$\psi(x,t+\Delta t) \approx e^{-iV(x)\Delta t/2} \, \mathcal{F}^{-1} \left[ e^{-i k^2 \Delta t/2} \, \mathcal{F}[\psi(x,t)] \right] \, e^{-iV(x)\Delta t/2}$$

This method is efficient because:

• the potential operator is simple in position space
• the kinetic operator is simple in momentum space
• FFT makes the transformations fast

---

Simulation Parameters

The main parameters used in the script are:

• spatial domain: $x \in [-10, 10]$
• number of grid points: $N = 1024$
• time step: $\Delta t = 0.005$
• total simulation time: $T = 20.0$

The spatial grid is defined uniformly, and the momentum grid is constructed using the FFT frequency convention.

These values are chosen to give stable and accurate results while keeping the code simple.

---

Quantities Computed

During the time evolution, the code computes several physically important quantities.

1. Probability Density

The probability density is

$$|\psi(x,t)|^2$$

This shows how the wavepacket is distributed in space at each time.

2. Expectation Value of Position

The expectation value of position is

$$\langle x(t) \rangle = \int x |\psi(x,t)|^2 \, dx$$

This quantity tracks the motion of the center of the wavepacket.

3. Normalization

The normalization of the wavefunction is

$$N(t) = \int |\psi(x,t)|^2 \, dx$$

For a physically correct and numerically stable simulation, this should remain close to 1.

4. Normalization Error

To check numerical accuracy, the code also computes

$$|N(t) - 1|$$

A very small value indicates that probability is well conserved during time evolution.

---

Expected Physical Behavior

For a Gaussian packet initially displaced from equilibrium in a harmonic potential, the packet oscillates back and forth in time.

The expectation value of position should show approximately sinusoidal behavior, consistent with harmonic motion. In other words, the center of the packet behaves in a way that reflects the classical oscillator, while the full evolution is still governed by quantum mechanics.

Because the harmonic oscillator is a smooth and well-behaved potential, it is a good system for testing the Split-Step Fourier method.

---

Output Plots

The script generates the following plots:

Initial Gaussian Wavepacket

Wavepacket Evolution Snapshots

Expectation Value of Position

Probability Conservation

Normalization Error

These plots help verify both the physical correctness of the simulation and the numerical stability of the method.

---

File Structure

The repository currently contains:

simulation.py
README.md
requirements.txt

• simulation.py contains the full implementation of the numerical method
• README.md explains the physics and the numerical approach
• requirements.txt lists the Python packages needed to run the code

---

How to Run

Install the required Python packages:

pip install -r requirements.txt

Then run the script:

python simulation.py

The code will generate the plots automatically.

---

Why This Project Matters

This project is useful as a simple computational physics implementation because it combines:

• quantum mechanics
• numerical simulation
• Fourier methods
• stability checking
• physical interpretation of results

It also provides a clear example of how numerical methods can be used to study the time evolution of quantum systems.

---

Possible Extensions

Some possible future improvements are:

• compare results for different time steps
• compare results for different spatial grid sizes
• compute the expectation value of momentum
• compute the energy expectation value
• compare the numerical solution with analytical harmonic oscillator results
• simulate coherent states
• apply the same method to other potentials such as barriers or double wells

---

Requirements

The code uses:

• numpy
• matplotlib

These are listed in requirements.txt.

---

Author

Umayr Utmaan

This repository was created as part of a computational physics study on numerical wavepacket dynamics in quantum systems.