Last Updated on June 10, 2025 by Sushanta Barman
Imaginary time evolution is a powerful computational technique widely used in quantum mechanics to find the ground state of a quantum system.
The idea stems from a mathematical transformation where time is treated as an imaginary quantity, allowing the evolution of a quantum state to “filter out” its ground state. This method is useful for systems where the ground state is not analytically obtainable, such as complex many-body systems or quantum fields.
In this article, we will explore the concept of imaginary time evolution, its mathematical formulation, and the steps to apply this technique to find the ground state of a given Hamiltonian.
Mathematical Background
Time-dependent Schrödinger equation
To understand imaginary time evolution, let’s start by recalling the time-dependent Schrödinger equation (TDSE),
\[i\hbar \frac{\partial \psi(x, t)}{\partial t} = \hat{H} \psi(x, t),\]
where, \(\psi(x, t)\) is the wave function of the quantum system that describes its state, \(\hbar\) is the reduced Planck’s constant, and \(\hat{H}\) is the Hamiltonian operator representing the total energy of the system.
The general solution to the Schrödinger equation can be written as a linear combination of the eigenstates \(\phi_n(x)\) of the Hamiltonian:
\[\psi(x, t) = \sum_{n} c_n \phi_n(x) e^{-i E_n t / \hbar},\]
where \(c_n\) are coefficients determined by the initial conditions of the system, and \(E_n\) are the eigenenergies corresponding to the eigenstates \(\phi_n(x)\).
Imaginary Time Evolution
To find the ground state, we consider a transformation from real time \(t\) to imaginary time \(\tau = it\). Substituting \(t = -i\tau\) into the TDSE, we obtain:
\[-\hbar \frac{\partial \psi(x, \tau)}{\partial \tau} = \hat{H} \psi(x, \tau),\]
or equivalently:
\[\frac{\partial \psi(x, \tau)}{\partial \tau} = -\frac{1}{\hbar} \hat{H} \psi(x, \tau).\]
This is known as the imaginary time Schrödinger equation. The solution to this equation takes the form:
\[\psi(x, \tau) = \sum_{n} c_n \phi_n(x) e^{-E_n \tau / \hbar}.\]
Ground State Filtering Mechanism
In the expression for \(\psi(x, \tau)\), each term involves an exponential factor \(e^{-E_n \tau / \hbar}\). As \(\tau \to \infty\), the term corresponding to the smallest eigenvalue \(E_0\) (the ground state energy) will dominate because it decays slower than all other terms. Hence, in the limit of large \(\tau\):
\[\psi(x, \tau) \approx c_0 \phi_0(x) e^{-E_0 \tau / \hbar}.\]
Therefore, by evolving a trial wave function \(\psi(x, 0)\) in imaginary time, it converges to the ground state \(\phi_0(x)\) up to a normalization factor.
Normalization and Repetition
To ensure numerical stability and prevent the wave function from becoming arbitrarily small, we normalize the wave function after each small step in \(\tau\). The normalized wave function at each step can be written as:
\[\psi(x, \tau + \Delta \tau) = \frac{\psi(x, \tau) – \frac{\Delta \tau}{\hbar} \hat{H} \psi(x, \tau)}{\|\psi(x, \tau) – \frac{\Delta \tau}{\hbar} \hat{H} \psi(x, \tau)\|}.\]
This process is repeated iteratively until convergence is achieved, where the state \(\psi(x, \tau)\) no longer changes significantly, indicating it has reached the ground state.
Algorithm for Imaginary Time Evolution
To implement the imaginary time evolution numerically, we use the following steps:
- Initialize the Wave Function:
Start with a trial wave function \(\psi(x, 0)\) that has non-zero overlap with the ground state \(\phi_0(x)\). This could be a random function or a physically motivated guess. - Discretize Time and Space:
Discretize the time evolution parameter \(\tau\) into small steps \(\Delta \tau\). Similarly, discretize the spatial domain \(x\) if working in a finite difference scheme. - Iterative Evolution:
Evolve the wave function in small imaginary time steps,
\[ \psi(x, \tau + \Delta \tau) = \frac{\psi(x, \tau) – \frac{\Delta \tau}{\hbar} \hat{H} \psi(x, \tau)}{\|\psi(x, \tau) – \frac{\Delta \tau}{\hbar} \hat{H} \psi(x, \tau)\|}. \] - Normalize the Wave Function:
Normalize \(\psi(x, \tau)\) after each step to avoid divergence or numerical overflow. - Check Convergence:
Continue the evolution until the change in \(\psi(x, \tau)\) between successive iterations is below a chosen threshold, indicating convergence to the ground state.
Example 1: Finding the Ground State of a 1D Quantum Harmonic Oscillator
In this section, the ground state of a one-dimensional quantum harmonic oscillator is determined by evolving a trial wave function in imaginary time, following the method outlined previously.
The system is described by the well-known harmonic potential \(V(x) = \frac{1}{2} m \omega^2 x^2\), and its corresponding Hamiltonian operator is given by \(\hat{H} = -\frac{\hbar^2}{2m} \frac{d^2}{dx^2} + V(x)\), which includes both kinetic and potential energy contributions.
Trail wave function
To begin the numerical process, a trial wave function is taken as:
\[\psi(x, 0) = e^{-x^2/2\sigma^2} \sin(kx)\].
This function is intentionally constructed to have non-zero overlap with both the ground state and some excited states of the system, enabling the filtering mechanism of imaginary time evolution to effectively isolate the lowest-energy state over repeated iterations.
As the wave function evolves in imaginary time, higher-energy components decay faster than the ground state component. Once convergence is achieved, the ground state energy is calculated using the expectation value of the Hamiltonian with respect to the evolved wave function:
\[E = \langle \psi | \hat{H} | \psi \rangle.\]
Analytical solution
Once the ground state wave function is obtained numerically, it is compared with the well-known analytical solution for the ground state, given by
\[\psi_0(x) = \left( \frac{m \omega}{\pi \hbar} \right)^{1/4} e^{-m \omega x^2 / (2 \hbar)},\]
and the corresponding ground state energy is given by the well-known expression,
\[E_0 = \frac{1}{2} \hbar \omega.\]
This provides a benchmark to validate the accuracy of the numerical solution obtained through imaginary time evolution.
| Parameter | Symbol | Value |
|---|---|---|
| Reduced Planck constant | \(\hbar\) | 1.0 |
| Mass of particle | \(m\) | 1.0 |
| Angular frequency | \(\omega\) | 1.0 |
| Number of spatial grid points | \(N_{\text{points}}\) | 500 |
| Spatial domain minimum | \(x_{\text{min}}\) | -7.0 |
| Spatial domain maximum | \(x_{\text{max}}\) | 7.0 |
| Width of Gaussian | \(\sigma\) | 2.0 |
| Wave vector | \(k\) | 5.0 |
| Imaginary time step | \(\Delta\tau\) | 0.0005 |
| Number of time steps | \(n_{\text{steps}}\) | 200,000 |
Python code
A Python code is provided below to compute the ground state energy and wave function starting from a trial wave function. All simulation parameters used in the code are listed in Table I. The simulation results are shown in Figure 1.
import numpy as np
import matplotlib.pyplot as plt
# Parameters
hbar = 1.0
m = 1.0
omega = 1.0
# 1. Spatial grid and harmonic oscillator potential
N_points = 500
x_min, x_max = -7.0, 7.0
x = np.linspace(x_min, x_max, N_points)
dx = x[1] - x[0]
V = 0.5 * m * omega**2 * x**2
# 2. Finite-difference Hamiltonian
main_diag = np.full(N_points, -2.0)
off_diag = np.full(N_points - 1, 1.0)
kinetic = -(hbar**2) / (2*m*dx**2) * (np.diag(main_diag) + np.diag(off_diag, 1) + np.diag(off_diag, -1))
potential = np.diag(V)
H = kinetic + potential
# 3. Gaussian-modulated sinusoidal initial wave function
sigma = 2.0
k = 5.0 # frequency of sinusoidal modulation
psi_init = np.exp(-x**2/(2*sigma**2)) * np.sin(k * x)
psi_init /= np.linalg.norm(psi_init)
psi = psi_init.copy()
# 4. Imaginary time evolution
dtau = 0.0005
n_steps = 200000
for i in range(n_steps):
psi -= dtau / hbar * H.dot(psi)
psi /= np.linalg.norm(psi)
# 5. Calculate ground state energy
E = np.dot(psi, H.dot(psi))
# Analytical ground state
psi_analytical = (m*omega/np.pi/hbar)**0.25 * np.exp(-m*omega*x**2/(2*hbar))
psi_analytical /= np.linalg.norm(psi_analytical)
E_analytical = 0.5 * hbar * omega
# Plotting
plt.figure(figsize=(8, 5))
plt.plot(x, psi_init, '-.g', label="Initial Wave Function", lw=2)
plt.plot(x, psi, '-b', label="Numerical Ground State", lw=2.5)
plt.plot(x, psi_analytical, '--r', label="Analytical Ground State", lw=2)
plt.xlabel(r"$x$" , fontsize=20)
plt.ylabel(r"$\psi(x)$", fontsize=20)
plt.title("Imaginary Time Evolution: 1D Quantum Harmonic Oscillator")
plt.legend(fontsize=13)
plt.tick_params(which="major", axis="both", direction="in", top=True, right=True, length=5, width=1, labelsize=16)
plt.ylim(-0.15,0.2)
#plt.grid()
plt.tight_layout()
plt.show()
print(f"Numerical ground state energy: {E:.6f}")
print(f"Analytical ground state energy: {E_analytical:.6f}")
About the Python code
The above Python code efficiently finds the ground state of a 1D quantum harmonic oscillator using imaginary time evolution. The finite-difference method accurately constructs the Hamiltonian, and the iterative normalization ensures numerical stability.
By choosing a well-shaped initial wave function, the algorithm rapidly filters out higher-energy states, converging to the ground state within a reasonable number of steps.
Key aspects include proper spatial discretization, normalization at each step, and the comparison with the analytical solution, which verifies the accuracy and effectiveness of the method.
Animation of the imaginary time evolution with iterations
Below is an animation video showing the evolution of the wave function and ground state energy with iterations during imaginary time evolution.
Here is the Python code for the above animation:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# Parameters
hbar = 1.0
m = 1.0
omega = 1.0
# 1. Spatial grid and harmonic oscillator potential
N_points = 500
x_min, x_max = -7.0, 7.0
x = np.linspace(x_min, x_max, N_points)
dx = x[1] - x[0]
V = 0.5 * m * omega**2 * x**2
# 2. Finite-difference Hamiltonian
main_diag = np.full(N_points, -2.0)
off_diag = np.full(N_points - 1, 1.0)
kinetic = -(hbar**2) / (2*m*dx**2) * (np.diag(main_diag) + np.diag(off_diag, 1) + np.diag(off_diag, -1))
potential = np.diag(V)
H = kinetic + potential
# 3. Gaussian-modulated sinusoidal initial wave function
sigma = 2.0
k = 5.0 # frequency of sinusoidal modulation
psi_init = np.exp(-x**2/(2*sigma**2)) * np.sin(k * x)
psi_init /= np.linalg.norm(psi_init)
psi = psi_init.copy()
# Analytical ground state
psi_analytical = (m*omega/np.pi/hbar)**0.25 * np.exp(-m*omega*x**2/(2*hbar))
psi_analytical /= np.linalg.norm(psi_analytical)
# Animation parameters
dtau = 0.0005
total_steps = 80000
frames = 400 # Total frames in the animation
steps_per_frame = total_steps // frames
# Store wave function evolution and energies for animation
wave_evolution = []
energies = []
psi_temp = psi_init.copy()
for i in range(frames):
for _ in range(steps_per_frame):
psi_temp -= dtau / hbar * H.dot(psi_temp)
psi_temp /= np.linalg.norm(psi_temp)
wave_evolution.append(psi_temp.copy())
energy = np.dot(psi_temp, H.dot(psi_temp))
energies.append(energy)
# Set up animation plot
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlabel(r"$x$", fontsize=20)
ax.set_ylabel(r"$\psi(x)$", fontsize=20)
# No static title here
ax.set_xlim(-7.0, 7.0)
ax.set_ylim(-0.15, 0.2)
ax.tick_params(which="major", axis="both", direction="in", top=True, right=True, length=5, width=1, labelsize=14)
line_num, = ax.plot([], [], '-b', lw=2.5, label="Numerical Wave Function")
line_init, = ax.plot([], [], '-.g', lw=2, label="Initial Wave Function")
line_analytical, = ax.plot([], [], '--r', lw=2, label="Analytical Ground State")
ax.legend(loc=1, fontsize=12)
# Add "MatterWaveX.Com" text at the lower right (location 4)
ax.text(0.98, 0.04, "MatterWaveX.Com", fontsize=13, color="gray",
ha='right', va='bottom', transform=ax.transAxes, alpha=0.7)
def init():
line_num.set_data([], [])
line_init.set_data(x, psi_init)
line_analytical.set_data(x, psi_analytical)
ax.set_title('') # Clear the title at init
return line_num, line_init, line_analytical
def animate(i):
line_num.set_data(x, wave_evolution[i])
iteration = (i + 1) * steps_per_frame
energy = energies[i]
ax.set_title(f"Iteration: {iteration} Ground State Energy: {energy:.6f}", fontsize=16)
return line_num, line_init, line_analytical
ani = animation.FuncAnimation(
fig, animate, init_func=init, frames=frames, interval=200, blit=True
)
plt.show()
# To save the animation as MP4:
ani.save("imaginary_time_evolution.mp4", writer="ffmpeg", fps=30)
Practical Applications
Imaginary time evolution is widely used in:
- Quantum Monte Carlo Simulations: To find ground states of complex many-body systems.
- Density Functional Theory (DFT): For minimizing the energy functional to find the ground state electron density.
- Condensed Matter Physics: To study properties of strongly correlated systems, like Bose-Einstein condensates.
Conclusion
Imaginary time evolution offers a straightforward and effective approach to finding the ground state of a quantum system, particularly in cases where direct analytical solutions are not possible.
By evolving a trial wave function in imaginary time, the technique systematically suppresses higher-energy components, isolating the ground state.
It is a critical tool in computational quantum mechanics and is fundamental for research in areas ranging from condensed matter physics to quantum chemistry.
By implementing the steps outlined in this article, you can use imaginary time evolution to find the ground state of any quantum system numerically, provided you have a suitable Hamiltonian and initial trial wave function.
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.
