Last Updated on June 15, 2025 by Sushanta Barman
Imaginary Time Evolution (ITE) is a computational technique extensively used in quantum mechanics to find the ground states of quantum systems numerically. It involves evolving an initial trial wave function in imaginary time, effectively damping higher energy components and isolating the ground state. The key idea relies on the mathematical transformation where the real-time Schrödinger equation is analytically continued into imaginary time, transforming the dynamics into a diffusion-like process.
The objective of this article is to provide a comprehensive understanding of how the ITE method can be generalized to compute excited states. This is particularly important because excited states offer valuable insights into quantum systems, their spectra, and transitions relevant in quantum technologies.
In this article, we discuss the fundamentals of ITE. The method for computing excited states using orthogonalization is presented. A practical implementation is demonstrated for a two-dimensional harmonic oscillator. The simulation results are interpreted and analyzed, and compared with the corresponding analytical solutions.
Theoritical background of Imaginary Time Evolution
Imaginary Time Evolution originates from the time-dependent Schrödinger equation (TDSE):
\[i\hbar \frac{\partial \psi(x, t)}{\partial t} = \hat{H}\psi(x, t)\]
where \(i\) is the imaginary unit, \(\hbar\) is the reduced Planck’s constant, \(\psi(x, t)\) is the wave function, and \(\hat{H}\) is the Hamiltonian operator of the system.
By substituting imaginary time \(\tau = it\), the equation becomes:
\[\frac{\partial \psi(x, \tau)}{\partial \tau} = -\frac{1}{\hbar}\hat{H}\psi(x, \tau)\]
where \(\tau\) is the imaginary time variable, and the equation now describes evolution in imaginary rather than real time.
This transformation converts the quantum dynamics into a diffusion equation form, where the wave function evolution dampens higher-energy components exponentially as follows:
\[\psi(x, \tau) = e^{-\hat{H}\tau/\hbar}\psi_{\text{trial}}(x)\]
where \(\psi_{\text{trial}}(x)\) is the initial trial wave function, and the exponential operator suppresses higher-energy eigenstates as \(\tau\) increases.
To compute excited states, the process involves iteratively orthogonalizing the wave function against previously obtained lower-energy states at each step, ensuring the convergence to higher-energy eigenstates.
Here, orthogonalization refers to the mathematical procedure of ensuring that the evolving wave function remains orthogonal (i.e., has zero overlap) with all lower-energy eigenstates previously computed, so that it can converge to the next higher-energy state.
Below, we demonstrate how this method can be employed to obtain both the ground and excited states of a two-dimensional quantum harmonic oscillator.
Computing Excited States of a Two-Dimensional Harmonic Oscillator
The two-dimensional harmonic oscillator is described by the Hamiltonian:
\[\hat{H} = -\frac{\hbar^2}{2m}\left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\right) + \frac{1}{2}m\omega^2(x^2 + y^2)\]
Analytical energy eigenvalues for the harmonic oscillator are known explicitly:
\[E_{n_x, n_y} = \hbar\omega(n_x + n_y + 1)\]
where \(n_x\) and \(n_y\) are quantum numbers.
Simulation Steps
Here we outline the principal steps used in the simulation for computing the ground and excited states. The approach is general and can be adapted to a wide class of Hamiltonians.
Discretization of Space
The simulation domain is discretized on a square grid:
- Define \(x\) and \(y\) coordinates over a uniform mesh.
- Construct the 2D meshgrid \(X, Y\) for spatial coordinates.
- The spatial step sizes \(\Delta x\) and \(\Delta y\) are computed from the grid. In our case, \(\Delta x=\Delta y\).
Definition of the Potential
The 2D harmonic oscillator potential is defined as:
\[ V(x, y) = \frac{1}{2}m\omega^2 (x^2 + y^2)\]
Initialization of Trial Wave Functions
Ground State: Start with a normalized Gaussian trial function,
\[\psi_0^{\text{trial}}(x, y) = \exp\left(-\frac{x^2 + y^2}{2\sigma_0^2}\right)\]
The distribution of \(\psi_0^{\text{trial}}(x, y)\) is shown in Figure 1.
Excited States: Use physically motivated functions (e.g., \(X\exp(-r^2/2\sigma_0^2)\) for \((1,0)\)), ensuring orthogonality to lower states.
Imaginary Time Evolution Loop
The evolution is governed by the imaginary time Schrödinger equation:
\[ \frac{\partial \psi}{\partial \tau} = -\frac{1}{\hbar}\hat{H}\psi,\]
where \(\tau = it\), and \(\hat{H}\) is the Hamiltonian.
At each time step:
1. Apply Hamiltonian:
- Compute the Laplacian (kinetic term) using a 5-point finite-difference stencil.
- Form the Hamiltonian action: \( H\psi = -\frac{\hbar^2}{2m}\nabla^2\psi + V(x, y)\psi\)
2. Update Wave Function: Use forward-Euler step as, \(\psi \leftarrow \psi – \frac{dtau}{\hbar}H\psi\)
3. Orthogonalization:
- For excited states, project out components along previously obtained lower-energy states: \( \psi \leftarrow \psi – \sum_k \langle \psi_k | \psi \rangle \psi_k\)
- Ensures convergence to a state orthogonal to all lower states.
4. Normalization: Normalize \(\psi\) at each iteration as,
\[\psi \leftarrow \frac{\psi}{\sqrt{\sum_{i,j} |\psi_{i,j}|^2 dx^2}}\]
Energy Calculation and Convergence Monitoring
- At each iteration, compute the expectation value of energy: \(E = \langle \psi | \hat{H} | \psi \rangle = \sum_{i,j} \psi^*_{i,j} \left[ -\frac{\hbar^2}{2m}\nabla^2\psi_{i,j} + V_{i,j}\psi_{i,j} \right] dx^2\)
- Track the convergence of the energy to ensure the state has stabilized.
Iterative Extraction of Excited States
- Ground state (\(n_x,n_y=0,0\)) is found with no orthogonalization.
- First excited state (\(n_x,n_y=1,0\)) is evolved with projection against the ground state.
- Second and higher excited states are obtained by projecting out all lower-lying converged states at each step.
Visualization and Analysis
After convergence:
- Plot the convergence history of energies for all computed states.
- Visualize the spatial structure (\(|\psi(x, y)|\)) of the trial and converged wave functions using 3D surface plots.
- Compare numerical energies with the analytic values: \(E_{n_x, n_y} = (n_x + n_y + 1)\hbar\omega\)
The above algorithm combines finite-difference discretization, imaginary-time propagation, Gram-Schmidt orthogonalization, and normalization to efficiently obtain ground and excited state wave functions of the 2D quantum harmonic oscillator. This approach is generalizable to many quantum systems where analytic solutions are not available.
Simulation Parameters
The numerical parameters used for our simulation are summarized in the table below:
| Parameter | Value |
|---|---|
| Grid size (N) | 200 |
| Spatial extent (x, y) | [-5, 5] |
| Mass (m) | 1.0 |
| Angular frequency (ω) | 1.0 |
| Reduced Planck’s constant (ℏ) | 1.0 |
| Time step (Δτ) | 0.001 |
| Number of iterations | 8000 |
| Gaussian width (σ) | 2.0 |
Python Code
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation
# =============================================================================
# 0. General Parameters
# =============================================================================
hbar = 1.0
m = 1.0
omega = 1.0
# Grid resolution
N = 200
x_min, x_max = -5.0, 5.0
y_min, y_max = x_min, x_max
x = np.linspace(x_min, x_max, N)
y = np.linspace(y_min, y_max, N)
dx = x[1] - x[0]
X, Y = np.meshgrid(x, y, indexing='ij')
# 2D harmonic oscillator potential: V(x,y) = (1/2) m ω^2 (x^2 + y^2)
V = 0.5 * m * omega**2 * (X**2 + Y**2)
# =============================================================================
# 1. Helper Functions
# =============================================================================
def normalize(psi, dx):
"""
Normalize a 2D wave function psi on a uniform square grid with spacing dx.
"""
norm = np.sqrt(np.sum(np.abs(psi)**2) * dx * dx)
return psi / norm
def laplacian_2d(psi, dx):
"""
5‐point finite‐difference Laplacian for a 2D array psi.
Dirichlet boundary assumed (psi=0 outside), so we only fill the interior.
"""
lap = np.zeros_like(psi, dtype=np.complex128)
lap[1:-1, 1:-1] = (
psi[2: , 1:-1] + psi[0:-2, 1:-1] +
psi[1:-1, 2: ] + psi[1:-1, 0:-2] -
4.0 * psi[1:-1, 1:-1]
) / dx**2
return lap
def project_out(psi, orthonormal_states, dx):
"""
Subtract from psi its components along each state in orthonormal_states.
Each state in orthonormal_states is assumed to be normalized.
Return the projected‐out psi.
"""
psi_proj = psi.copy()
for phi in orthonormal_states:
overlap = np.sum(np.conj(phi) * psi_proj) * dx * dx
psi_proj = psi_proj - overlap * phi
return psi_proj
# =============================================================================
# 2. Analytical Ground State (for comparison)
# =============================================================================
psi_analytical = (m * omega / (np.pi * hbar))**0.5 * np.exp(-0.5 * m * omega * (X**2 + Y**2) / hbar)
psi_analytical = normalize(psi_analytical, dx)
# =============================================================================
# 3. Imaginary‐Time Evolution Function
# =============================================================================
def imaginary_time_evolve(
psi_initial,
V,
dx,
dtau,
total_steps,
orthonormal_states=None
):
"""
Perform imaginary-time evolution on psi_initial to find the lowest‐energy
state orthogonal to any in 'orthonormal_states'.
"""
psi = psi_initial.copy()
if orthonormal_states is None:
orthonormal_states = []
E_history = []
kinetic_prefactor = - (hbar**2) / (2.0 * m)
for step in range(1, total_steps + 1):
lap_psi = laplacian_2d(psi, dx)
H_psi = kinetic_prefactor * lap_psi + V * psi
psi = psi - (dtau / hbar) * H_psi
if orthonormal_states:
psi = project_out(psi, orthonormal_states, dx)
psi = normalize(psi, dx)
E_loc = np.sum(np.conj(psi) * (kinetic_prefactor * laplacian_2d(psi, dx) + V * psi)) * dx * dx
E_history.append(np.real(E_loc))
return psi, E_history
# =============================================================================
# 4. Set Up & Find ψ₀ (Ground State: (n_x, n_y) = (0,0))
# =============================================================================
sigma0 = 2.0
psi0_trial = np.exp(-((X**2 + Y**2) / (2.0 * sigma0**2)))
psi0_trial = normalize(psi0_trial, dx)
dtau = 0.001
total_steps = 8000
print("Evolving to find the GROUND STATE (n_x, n_y) = (0, 0)...")
psi0, E0_hist = imaginary_time_evolve(
psi_initial = psi0_trial,
V = V,
dx = dx,
dtau = dtau,
total_steps = total_steps,
orthonormal_states = None
)
E0_numerical = E0_hist[-1]
print(f"Converged Ground‐State Energy (numerical): {E0_numerical:.6f}")
print(f"Analytic Ground‐State Energy: {1.0 * hbar * omega:.6f}\n")
# =============================================================================
# 5. Find ψ₁ (First Excited State: (n_x, n_y) = (1,0))
# =============================================================================
psi1_trial = X * np.exp(-((X**2 + Y**2) / (2.0 * sigma0**2)))
psi1_trial = psi1_trial - (np.sum(np.conj(psi0) * psi1_trial) * dx * dx) * psi0
psi1_trial = normalize(psi1_trial, dx)
print("Evolving to find the FIRST EXCITED STATE (n_x, n_y) = (1, 0)...")
psi1, E1_hist = imaginary_time_evolve(
psi_initial = psi1_trial,
V = V,
dx = dx,
dtau = dtau,
total_steps = total_steps,
orthonormal_states = [psi0]
)
E1_numerical = E1_hist[-1]
print(f"Converged First Excited‐State Energy (numerical): {E1_numerical:.6f}")
print(f"Analytic First Excited‐State Energy: {2.0 * hbar * omega:.6f}\n")
# =============================================================================
# 6. Find ψ₂ (Second Excited State: (n_x, n_y) = (0,1))
# =============================================================================
psi2_trial = Y * np.exp(-((X**2 + Y**2) / (2.0 * sigma0**2)))
psi2_trial = psi2_trial - (np.sum(np.conj(psi0) * psi2_trial) * dx * dx) * psi0
psi2_trial = psi2_trial - (np.sum(np.conj(psi1) * psi2_trial) * dx * dx) * psi1
psi2_trial = normalize(psi2_trial, dx)
print("Evolving to find the SECOND EXCITED STATE (n_x, n_y) = (0, 1)...")
psi2, E2_hist = imaginary_time_evolve(
psi_initial = psi2_trial,
V = V,
dx = dx,
dtau = dtau,
total_steps = total_steps,
orthonormal_states = [psi0, psi1]
)
E2_numerical = E2_hist[-1]
print(f"Converged Second Excited‐State Energy (numerical): {E2_numerical:.6f}")
print(f"Analytic Second Excited‐State Energy: {2.0 * hbar * omega:.6f}\n")
# =============================================================================
# 7. Find ψ₃ (Third Excited State: (n_x, n_y) = (2,0))
# =============================================================================
psi3_trial = (X**2) * np.exp(-((X**2 + Y**2) / (2.0 * sigma0**2)))
psi3_trial = psi3_trial - (np.sum(np.conj(psi0) * psi3_trial) * dx * dx) * psi0
psi3_trial = psi3_trial - (np.sum(np.conj(psi1) * psi3_trial) * dx * dx) * psi1
psi3_trial = psi3_trial - (np.sum(np.conj(psi2) * psi3_trial) * dx * dx) * psi2
psi3_trial = normalize(psi3_trial, dx)
print("Evolving to find the THIRD EXCITED STATE (n_x, n_y) = (2, 0)...")
psi3, E3_hist = imaginary_time_evolve(
psi_initial = psi3_trial,
V = V,
dx = dx,
dtau = dtau,
total_steps = total_steps,
orthonormal_states = [psi0, psi1, psi2]
)
E3_numerical = E3_hist[-1]
print(f"Converged Third Excited‐State Energy (numerical): {E3_numerical:.6f}\n")
# =============================================================================
# 8. Find ψ₄ (Fourth Excited State: (n_x, n_y) = (1,1))
# =============================================================================
psi4_trial = (X * Y) * np.exp(-((X**2 + Y**2) / (2.0 * sigma0**2)))
psi4_trial = psi4_trial - (np.sum(np.conj(psi0) * psi4_trial) * dx * dx) * psi0
psi4_trial = psi4_trial - (np.sum(np.conj(psi1) * psi4_trial) * dx * dx) * psi1
psi4_trial = psi4_trial - (np.sum(np.conj(psi2) * psi4_trial) * dx * dx) * psi2
psi4_trial = psi4_trial - (np.sum(np.conj(psi3) * psi4_trial) * dx * dx) * psi3
psi4_trial = normalize(psi4_trial, dx)
print("Evolving to find the FOURTH EXCITED STATE (n_x, n_y) = (1, 1)...")
psi4, E4_hist = imaginary_time_evolve(
psi_initial = psi4_trial,
V = V,
dx = dx,
dtau = dtau,
total_steps = total_steps,
orthonormal_states = [psi0, psi1, psi2, psi3]
)
E4_numerical = E4_hist[-1]
print(f"Converged Fourth Excited‐State Energy (numerical): {E4_numerical:.6f}\n")
# =============================================================================
# 9. Find ψ₅ (Fifth Excited State: (n_x, n_y) = (0,2))
# =============================================================================
psi5_trial = (Y**2) * np.exp(-((X**2 + Y**2) / (2.0 * sigma0**2)))
psi5_trial = psi5_trial - (np.sum(np.conj(psi0) * psi5_trial) * dx * dx) * psi0
psi5_trial = psi5_trial - (np.sum(np.conj(psi1) * psi5_trial) * dx * dx) * psi1
psi5_trial = psi5_trial - (np.sum(np.conj(psi2) * psi5_trial) * dx * dx) * psi2
psi5_trial = psi5_trial - (np.sum(np.conj(psi3) * psi5_trial) * dx * dx) * psi3
psi5_trial = psi5_trial - (np.sum(np.conj(psi4) * psi5_trial) * dx * dx) * psi4
psi5_trial = normalize(psi5_trial, dx)
print("Evolving to find the FIFTH EXCITED STATE (n_x, n_y) = (0, 2)...")
psi5, E5_hist = imaginary_time_evolve(
psi_initial = psi5_trial,
V = V,
dx = dx,
dtau = dtau,
total_steps = total_steps,
orthonormal_states = [psi0, psi1, psi2, psi3, psi4]
)
E5_numerical = E5_hist[-1]
print(f"Converged Fifth Excited‐State Energy (numerical): {E5_numerical:.6f}\n")
# =============================================================================
# 10. Convergence Histories for All States
# =============================================================================
plt.figure(figsize=(8, 6))
plt.plot(E0_hist, label=r"$\psi_{0}$ $(n_x,n_y)=(0,0)$")
plt.plot(E1_hist, label=r"$\psi_{1}$ $(n_x,n_y)=(1,0)$")
plt.plot(E2_hist, label=r"$\psi_{2}$ $(n_x,n_y)=(0,1)$")
plt.plot(E3_hist, label=r"$\psi_{3}$ $(n_x,n_y)=(2,0)$")
plt.plot(E4_hist, label=r"$\psi_{4}$ $(n_x,n_y)=(1,1)$")
plt.plot(E5_hist, label=r"$\psi_{5}$ $(n_x,n_y)=(0,2)$")
plt.hlines(1.0, 0, total_steps, colors='k', linestyles='--', label=r"$E=(0+0+1)\hbar\omega$")
plt.hlines(2.0, 0, total_steps, colors='r', linestyles='--', label=r"$E=(1+0+1)\hbar\omega = (0+1+1)\hbar\omega$")
plt.hlines(3.0, 0, total_steps, colors='b', linestyles='--', label=r"$E=(2+0+1)\hbar\omega = (1+1+1)\hbar\omega = (0+2+1)\hbar\omega$")
plt.xlabel("Iteration", fontsize=14)
plt.ylabel("Energy", fontsize=14)
plt.title("Convergence of Imaginary‐Time Evolution for First Six States")
plt.legend(fontsize=10)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.tight_layout()
plt.show()
# =============================================================================
# 11. 3×3 “Snapshot” Plot Including ψ₀_trial and ψ₀ … ψ₅ with (n_x, n_y)
# =============================================================================
fig = plt.figure(figsize=(15, 12))
# Row 1, Col 1: |ψ₀_trial| (trial for ground state)
ax1 = fig.add_subplot(3, 3, 1, projection='3d')
Z1 = np.abs(psi0_trial)
surf1 = ax1.plot_surface(
X, Y, Z1,
cmap='viridis', linewidth=0, antialiased=False, edgecolor='none'
)
ax1.set_title("Trial for Ground State $|\psi_{0,\mathrm{trial}}|$ $(n_x,n_y)=(0,0)$", fontsize=12)
ax1.set_xlabel("x", fontsize=12)
ax1.set_ylabel("y", fontsize=12)
ax1.set_zlabel(r"$|\psi_{\!0,\mathrm{trial}}|$", fontsize=12)
ax1.set_zlim(0, np.max(Z1) * 1.1)
ax1.tick_params(axis='x', labelsize=10)
ax1.tick_params(axis='y', labelsize=10)
ax1.tick_params(axis='z', labelsize=10)
# Row 1, Col 2: |ψ₀| (Converged Ground State)
ax2 = fig.add_subplot(3, 3, 2, projection='3d')
Z2 = np.abs(psi0)
surf2 = ax2.plot_surface(
X, Y, Z2,
cmap='viridis', linewidth=0, antialiased=False, edgecolor='none'
)
ax2.set_title("Converged Ground State $|\psi_{0}|$ $(n_x,n_y)=(0,0)$", fontsize=12)
ax2.set_xlabel("x", fontsize=12)
ax2.set_ylabel("y", fontsize=12)
ax2.set_zlabel(r"$|\psi_{0}|$", fontsize=12)
ax2.set_zlim(0, np.max(np.abs(psi_analytical)) * 1.1)
ax2.tick_params(axis='x', labelsize=10)
ax2.tick_params(axis='y', labelsize=10)
ax2.tick_params(axis='z', labelsize=10)
# Row 1, Col 3: |ψ₁| (First Excited)
ax3 = fig.add_subplot(3, 3, 3, projection='3d')
Z3 = np.abs(psi1)
surf3 = ax3.plot_surface(
X, Y, Z3,
cmap='viridis', linewidth=0, antialiased=False, edgecolor='none'
)
ax3.set_title("1st Excited State $|\psi_{1}|$ $(n_x,n_y)=(1,0)$", fontsize=12)
ax3.set_xlabel("x", fontsize=12)
ax3.set_ylabel("y", fontsize=12)
ax3.set_zlabel(r"$|\psi_{1}|$", fontsize=12)
ax3.set_zlim(0, np.max(np.abs(psi_analytical)) * 1.1)
ax3.tick_params(axis='x', labelsize=10)
ax3.tick_params(axis='y', labelsize=10)
ax3.tick_params(axis='z', labelsize=10)
# Row 2, Col 1: |ψ₂| (Second Excited)
ax4 = fig.add_subplot(3, 3, 4, projection='3d')
Z4 = np.abs(psi2)
surf4 = ax4.plot_surface(
X, Y, Z4,
cmap='viridis', linewidth=0, antialiased=False, edgecolor='none'
)
ax4.set_title("2nd Excited State $|\psi_{2}|$ $(n_x,n_y)=(0,1)$", fontsize=12)
ax4.set_xlabel("x", fontsize=12)
ax4.set_ylabel("y", fontsize=12)
ax4.set_zlabel(r"$|\psi_{2}|$", fontsize=12)
ax4.set_zlim(0, np.max(np.abs(psi_analytical)) * 1.1)
ax4.tick_params(axis='x', labelsize=10)
ax4.tick_params(axis='y', labelsize=10)
ax4.tick_params(axis='z', labelsize=10)
# Row 2, Col 2: |ψ₃| (Third Excited)
ax5 = fig.add_subplot(3, 3, 5, projection='3d')
Z5 = np.abs(psi3)
surf5 = ax5.plot_surface(
X, Y, Z5,
cmap='viridis', linewidth=0, antialiased=False, edgecolor='none'
)
ax5.set_title("3rd Excited State $|\psi_{3}|$ $(n_x,n_y)=(2,0)$", fontsize=12)
ax5.set_xlabel("x", fontsize=12)
ax5.set_ylabel("y", fontsize=12)
ax5.set_zlabel(r"$|\psi_{3}|$", fontsize=12)
ax5.set_zlim(0, np.max(np.abs(psi_analytical)) * 1.1)
ax5.tick_params(axis='x', labelsize=10)
ax5.tick_params(axis='y', labelsize=10)
ax5.tick_params(axis='z', labelsize=10)
# Row 2, Col 3: |ψ₄| (Fourth Excited)
ax6 = fig.add_subplot(3, 3, 6, projection='3d')
Z6 = np.abs(psi4)
surf6 = ax6.plot_surface(
X, Y, Z6,
cmap='viridis', linewidth=0, antialiased=False, edgecolor='none'
)
ax6.set_title("4th Excited State $|\psi_{4}|$ $(n_x,n_y)=(1,1)$", fontsize=12)
ax6.set_xlabel("x", fontsize=12)
ax6.set_ylabel("y", fontsize=12)
ax6.set_zlabel(r"$|\psi_{4}|$", fontsize=12)
ax6.set_zlim(0, np.max(np.abs(psi_analytical)) * 1.1)
ax6.tick_params(axis='x', labelsize=10)
ax6.tick_params(axis='y', labelsize=10)
ax6.tick_params(axis='z', labelsize=10)
# Row 3, Col 1: |ψ₅| (Fifth Excited)
ax7 = fig.add_subplot(3, 3, 7, projection='3d')
Z7 = np.abs(psi5)
surf7 = ax7.plot_surface(
X, Y, Z7,
cmap='viridis', linewidth=0, antialiased=False, edgecolor='none'
)
ax7.set_title("5th Excited State $|\psi_{5}|$ $(n_x,n_y)=(0,2)$", fontsize=12)
ax7.set_xlabel("x", fontsize=12)
ax7.set_ylabel("y", fontsize=12)
ax7.set_zlabel(r"$|\psi_{5}|$", fontsize=12)
ax7.set_zlim(0, np.max(np.abs(psi_analytical)) * 1.1)
ax7.tick_params(axis='x', labelsize=10)
ax7.tick_params(axis='y', labelsize=10)
ax7.tick_params(axis='z', labelsize=10)
# Keep the last two subplots blank
ax8 = fig.add_subplot(3, 3, 8)
ax8.axis('off')
ax9 = fig.add_subplot(3, 3, 9)
ax9.axis('off')
plt.tight_layout()
plt.show()
# =============================================================================
# 12. Summary of Results
# =============================================================================
print("Summary of Converged Energies:")
print(f" Ground State ψ₀ (n_x,n_y)=(0,0): E₀ (numerical) = {E0_numerical:.6f}, E₀ (analytic) = {1.0*hbar*omega:.6f}")
print(f" 1st Excited State ψ₁ (n_x,n_y)=(1,0): E₁ (numerical) = {E1_numerical:.6f}, E₁ (analytic) = {2.0*hbar*omega:.6f}")
print(f" 2nd Excited State ψ₂ (n_x,n_y)=(0,1): E₂ (numerical) = {E2_numerical:.6f}, E₂ (analytic) = {2.0*hbar*omega:.6f}")
print(f" 3rd Excited State ψ₃ (n_x,n_y)=(2,0): E₃ (numerical) = {E3_numerical:.6f}")
print(f" 4th Excited State ψ₄ (n_x,n_y)=(1,1): E₄ (numerical) = {E4_numerical:.6f}")
print(f" 5th Excited State ψ₅ (n_x,n_y)=(0,2): E₅ (numerical) = {E5_numerical:.6f}")
Get the above code in MATLAB:
Simulation Outputs and Analysis
The above Python code computes six distinct wave functions and their corresponding energies for the two-dimensional quantum harmonic oscillator. These include the ground state and five excited states, each characterized by quantum numbers \((n_x, n_y)\):
- The ground state \((0,0)\) and the first excited state \((1,0)\) are presented in Figure 2.
- The second \((0,1)\) and third excited states \((2,0)\) are shown in Figure 3.
- The fourth \((1,1)\) and fifth excited states \((0,2)\) are illustrated in Figure 4.
Each state clearly demonstrates the characteristic nodal structures associated with the respective quantum numbers of the two-dimensional harmonic oscillator eigenstates.
Analysis of Simulation Results
- Ground State: Numerically converged smoothly to a Gaussian shape, exactly matching analytical predictions.
- Excited States: Wave functions displayed distinct nodal patterns (as observed from the plots in Figures 2–4), accurately reflecting the theoretical predictions for quantum numbers \((n_x, n_y)\).
Comparison with Analytical Solutions
The convergence of the energies for all computed states at different iteration steps, obtained via imaginary time evolution, is shown in Figure 5. It can be observed that the numerical energies are in good agreement with the analytical solutions.
- Ground state (0,0): Numerically computed value is very close to the exact analytical value of \(E_0 = 1 \hbar \omega\).
- First and second excited states (1,0 and 0,1): Matched analytical energy \(E= 2 \hbar \omega\).
- Higher excited states: Also aligned closely with analytical predictions, showing minimal numerical errors and high accuracy.
The convergence plots shown in Figure 5 clearly illustrate how the energy of each state stabilizes after sufficient iterations, indicating the reliability of the method.
Summary
Imaginary Time Evolution (ITE) provides an effective computational approach not only for ground states but also for systematically obtaining excited states through orthogonalization.
This technique, demonstrated using a 2D harmonic oscillator, yields results closely matching analytical solutions, validating its accuracy and applicability.
The robustness of this method opens up numerous possibilities in exploring quantum spectra, energy transitions, and intricate quantum behaviors across various complex systems relevant to advanced quantum mechanics research and technology.
Further Reading
[1] Sushanta Barman, “Finding Ground State of 2D Quantum Harmonic Oscillator Using Imaginary Time Evolution,” MatterWaveX.com, (accessed June 1, 2025).
[2] Sushanta Barman, “Imaginary Time Evolution vs Real Time Evolution,” MatterWaveX.com, (accessed June 1, 2025).
[3] Sushanta Barman, “Imaginary Time Evolution Method for Finding the Ground State,” MatterWaveX.com, (accessed June 1, 2025).
[4] Q.-X. Xie, Y. Song, Y. Liu, and Y. Yao, “Improved algorithms of quantum imaginary time evolution for ground-state preparation,” J. Chem. Theory Comput. 19, 503 (2023).
[5] A. K. Roy, “Ground and excited states of spherically symmetric potentials through an imaginary-time evolution method: Application to spiked harmonic oscillators,” J. Math. Chem. 52, 2645 (2014).
[6] J. Yang and T. I. Lakoba, “Accelerated imaginary-time evolution methods for the computation of solitary waves,” Stud. Appl. Math. 120, 265 (2008).
[7] M. E. Stroeks, J. Helsen, and B. M. Terhal, “Spectral estimation for Hamiltonians: A comparison between classical imaginary-time evolution and quantum real-time evolution,” New J. Phys. 24, 103024 (2022).
[8] M. Motta et al., “Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution,” Nat. Phys. 16, 205 (2020).
I am a senior research scholar in the Department of Physics at IIT Kanpur. My work focuses on ion-beam optics and matter-wave phenomena. I am also interested in emerging matter-wave technologies.
