Particle in a Box: Infinite Square Well

Learning Objectives

By the end of this section, you will be able to:

1. Set up the time-independent Schrödinger equation inside a hard-walled potential and read off the boundary conditions straight from the physics.
2. Solve the equation by hand and obtain the eigenfunctions $\psi_n(x) = \sqrt{2/L}\,\sin(n\pi x/L)$ and energies $E_n = n^2\pi^2\hbar^2/(2mL^2)$.
3. Explain why the energy spectrum is discrete and why the ground state cannot have $E = 0$.
4. Interpret the nodes of $\psi_n$ and use orthogonality $\int_0^L \psi_m\psi_n\,dx = \delta_{mn}$ to expand any state.
5. Compute a real example (electron in a 1-nm box) end to end, in eV and in photon-emission wavelengths.
6. Watch a superposition oscillate in time and identify the beat frequency $\omega_{mn} = (E_m-E_n)/\hbar$.
7. Verify the analytic result numerically by discretizing the Hamiltonian and diagonalizing it in Python and PyTorch.

The Big Picture

"The particle in a box is to quantum mechanics what the simple harmonic oscillator is to classical mechanics: the first place where everything works out cleanly and you can see all the moving parts."

Imagine a bead trapped between two impassable walls a distance $L$ apart, free to slide between them but unable to escape. Classically there is nothing to talk about — the bead has whatever energy you gave it and bounces forever. Quantum mechanically, this is the simplest non-trivial problem in physics, and it already contains three of the deepest surprises of the theory:

1. Energy is quantized. Only certain discrete energies are allowed: $E_1, 4E_1, 9E_1, 16E_1, \ldots$
2. The lowest energy is not zero. Even at "rest" the particle has zero-point energy $E_1 = \pi^2\hbar^2/(2mL^2)$.
3. Standing waves replace trajectories. The state is a shape, not a position. The shape has nodes — places where the particle is never found.

Setting Up the Box

The infinite square well is the potential

The infinite barrier outside $[0,L]$ is a mathematical idealization of a physical wall the particle can absolutely not penetrate. The wave function is forced to be zero outside the box, and by continuity it must also vanish on the boundary:

Inside the box $V = 0$, so the time-independent Schrödinger equation $\hat H\psi = E\psi$ reduces to:

That is a second-order linear ODE with constant coefficients — the kind of object you have been solving since first-year calculus. The only quantum part is what we will do with it at the boundary.

Solving the Schrödinger Equation

Rearrange the ODE into the familiar harmonic-oscillator form:

For any $E > 0$ we have $k > 0$ and the general solution is the familiar pair of sines and cosines:

Step 1 — Apply $\psi(0) = 0$

Plug in $x = 0$: $\psi(0) = A\sin(0) + B\cos(0) = B$. For this to vanish we need $B = 0$. The cosine piece is killed by the left wall.

Step 2 — Apply $\psi(L) = 0$

What is left is $\psi(x) = A\sin(kx)$. The right wall demands $A\sin(kL) = 0$. Setting $A = 0$ would zero out the wave function everywhere — that is the trivial solution and corresponds to no particle at all. The only useful way to satisfy the equation is to force the sine itself to vanish:

Negative $n$ just flips the overall sign of $\psi$ and doesn't give a new physical state. $n = 0$ gives $\psi \equiv 0$, so we throw it out. The allowed values of $k$ are therefore:

Boundary Conditions Force Quantization

Look at what just happened. The ODE itself doesn't care about $E$ — any positive energy gives a perfectly valid sine wave. It was the boundary conditions that selected a discrete set. Translating $k_n = n\pi/L$ back into energy via $k^2 = 2mE/\hbar^2$:

This is energy quantization, derived from nothing more sophisticated than "the wave has to fit between two walls." The intuition is exactly the same as for a guitar string: only certain wavelengths fit, and each allowed wavelength corresponds to a definite frequency — and through the de Broglie relation $p = \hbar k$, a definite momentum and energy.

Guitar string analogy. A string fixed at both ends vibrates only at frequencies whose half-wavelength fits an integer number of times into its length: $L = n\lambda/2$. The quantum particle in a box is the exact same equation, with the energy taking the role of the squared frequency. Quantum mechanics didn't invent discrete spectra — every musical instrument has one. The difference is that for the bead in the box, the "string" is the wave function itself.

Normalization and Orthogonality

The wave function still has an unknown amplitude $A$. We fix it by demanding total probability one:

With $\psi_n(x) = A\sin(n\pi x/L)$ the integral evaluates to:

which yields $A = \sqrt{2/L}$. So our normalized eigenstates are:

Orthogonality

Different eigenstates are orthogonal: for $m \neq n$,

Combined with normalization this is $\int_0^L \psi_m\psi_n\,dx = \delta_{mn}$. A quick way to see it: use the product-to-sum identity $\sin a\sin b = \tfrac{1}{2}[\cos(a-b)-\cos(a+b)]$. Both cosine integrals over $[0,L]$ vanish when $m \neq n$, because the cosine completes whole periods on the interval.

Interactive: Eigenstate Explorer

Move the sliders. n picks the quantum number; L stretches or shrinks the box. Watch what each change does:

• Increasing n adds nodes (interior zeros where $\psi = 0$) and pushes the energy up like $n^2$.
• Increasing L stretches the wave horizontally and lowers every energy as $1/L^2$.
• For even $n$, the box centre $x = L/2$ is a node — the particle is never found exactly in the middle. For odd $n$ it is an antinode (the most likely spot).

The Energy Ladder

The energies form a ladder whose rungs spread out quadratically:

Two things are worth pausing on. First, the gap $E_{n+1}-E_n = (2n+1)E_1$ grows linearly with $n$ — the higher you climb, the more energy each new rung costs. Second, therelative gap $(E_{n+1}-E_n)/E_n \to 2/n$ shrinks. In the limit of huge $n$ the spectrum looks effectively continuous — this is the correspondence principle: classical physics emerges smoothly from quantum mechanics at high quantum numbers.

Properties of the Eigenstates

Nodes

$\psi_n(x) = \sqrt{2/L}\sin(n\pi x/L)$ vanishes at $x = kL/n$ for $k = 0, 1, 2, \ldots, n$. The endpoints $x=0$ and $x=L$ are forced zeros (the walls), leaving n − 1 interior nodes. The number of nodes is the quickest fingerprint of which eigenstate you are looking at.

Parity around the centre

Define $u = x - L/2$. Then $\sin(n\pi x/L) = \sin(n\pi/2)\cos(n\pi u/L) + \cos(n\pi/2)\sin(n\pi u/L)$. For odd $n$ only the cosine term survives — the eigenstate is symmetric about $L/2$. For even $n$ only the sine survives — it is antisymmetric. This alternation is why even-n states have a node right at the centre.

Expectation values you can compute by hand

Worked Example: Electron in a 1 nm Box

Numbers grounded a model. Take an electron confined to a one-nanometre box — a reasonable cartoon of a quantum dot or a single benzene-ring-size cavity.

Superposition and Time Evolution

A single eigenstate evolves only in phase: $\Psi_n(x,t) = \psi_n(x)\,e^{-iE_n t/\hbar}$. The complex phase rotates but $|\Psi_n(x,t)|^2 = |\psi_n(x)|^2$ doesn't change. Stationary states are truly stationary — their probability density is frozen.

The interesting dynamics show up the moment you combine eigenstates with different energies. Suppose:

Square it. The cross term picks up $\cos\!\big((E_2-E_1)t/\hbar\big)$:

The first two pieces are static. The last piece oscillates at the beat (Bohr) frequency $\omega_{21} = (E_2-E_1)/\hbar$. Because $\psi_1$ is symmetric and $\psi_2$ is antisymmetric about $L/2$, the product $\psi_1\psi_2$ is antisymmetric: it is positive on one half of the box and negative on the other. So as the cosine oscillates, the probability density sloshes from one side to the other.

Numerical Solution in Python

It is a one-line story: discretize the Hamiltonian on a grid, diagonalize the matrix, read off the spectrum. The finite-difference rule turns the operator $-(\hbar^2/2m)\,d^2/dx^2$ into a tridiagonal matrix. Diagonalizing it returns approximate eigenvalues — the discrete $E_n$ — and discrete samples of the eigenfunctions. The boundary conditions $\psi(0)=\psi(L)=0$ are built in just by the choice of grid (we omit the wall points from the unknowns).

np.array([1.0, 2.0, 3.0])

n=1: numerical=9.8694, exact=9.8696

dx = 1.0 / 201 ≈ 0.004975

x = [0.005, 0.010, 0.015, ..., 0.995]

ψ''(x_5) ≈ (ψ_6 − 2ψ_5 + ψ_4) / dx²

H = [[1/dx², -1/2dx², 0, ...], [-1/2dx², 1/dx², -1/2dx², ...], ...]

energies[0] = lowest eigenvalue ≈ π²/(2L²) ≈ 4.9348

states[:, 0] = the ground-state wavefunction sampled at every x_i

1import numpy as np
2import matplotlib.pyplot as plt
3
4# -------- 1.  Physical setup ---------------------------------------------
5# Natural units: hbar = 1, m = 1.  Box length L is the only scale.
6L = 1.0
7N = 200              # number of interior grid points
8dx = L / (N + 1)     # spacing between adjacent grid points
9
10# Interior position grid: x_1, x_2, ..., x_N
11# The walls live at x_0 = 0 and x_{N+1} = L; they are NOT in the array.
12x = np.linspace(dx, L - dx, N)
13
14# -------- 2.  Build the Hamiltonian H ------------------------------------
15# Inside the box V(x) = 0, so H = -(1/2) d^2/dx^2.
16# Three-point finite difference for the second derivative:
17#   psi''(x_i) ~= ( psi_{i+1} - 2*psi_i + psi_{i-1} ) / dx^2
18# Drop psi_0 and psi_{N+1}: that *automatically* enforces psi(0)=psi(L)=0.
19main_diag = (1.0 / dx**2) * np.ones(N)              # = -(-2)/(2 dx^2) = 1/dx^2
20off_diag  = (-0.5 / dx**2) * np.ones(N - 1)         # = -(1)/(2 dx^2)
21H = (
22    np.diag(main_diag)
23    + np.diag(off_diag, k=+1)
24    + np.diag(off_diag, k=-1)
25)
26
27# -------- 3.  Diagonalize ------------------------------------------------
28# H is real symmetric -> np.linalg.eigh returns (eigvals, eigvecs)
29# sorted by increasing eigenvalue.  Each column eigvecs[:, n] is psi_n.
30energies, states = np.linalg.eigh(H)
31
32# -------- 4.  Compare to the analytic result E_n = n^2 pi^2 / (2 L^2) ----
33print("First five energy levels (numerical vs exact):")
34for n in range(1, 6):
35    exact = (n * np.pi) ** 2 / (2 * L**2)
36    err   = abs(energies[n - 1] - exact)
37    print(f"  n={n}:  numerical={energies[n-1]:8.4f}   "
38          f"exact={exact:8.4f}   error={err:.2e}")
39
40# -------- 5.  Normalize and plot the first three eigenstates -------------
41fig, ax = plt.subplots(figsize=(8, 4))
42for n in range(3):
43    psi = states[:, n]
44    # Discrete normalization: sum |psi_i|^2 * dx = 1
45    psi = psi / np.sqrt(np.sum(psi**2) * dx)
46    # Fix sign so psi_n(L/2) > 0 for odd n -- pure cosmetics
47    if psi[N // 2] < 0:
48        psi = -psi
49    ax.plot(x, psi, label=f"n={n+1}")
50
51ax.axhline(0, color="gray", lw=0.5)
52ax.set_xlabel("x")
53ax.set_ylabel(r"$\psi_n(x)$")
54ax.legend()
55plt.show()

The Same Idea in PyTorch

Why bother rewriting in PyTorch? Three good reasons. First, the operation is GPU-friendly: $\hat H$ as a matrix is exactly what GPUs are designed for. Second, PyTorch's $\texttt{linalg.eigh}$ works on batches, so you can solve the eigenproblem for many different potentials or many different box lengths in one call. Third, the same Hamiltonian shows up as a building block in physics-informed neural networks and neural quantum states — having the operator in PyTorch lets you wire it directly into a learning loop.

device = cuda  →  4× speedup when sweeping 1024 box lengths

energies, states = torch.linalg.eigh(H)  # H must be symmetric

energies[0].item() → 4.9347

Ls = [0.50, 0.60, 0.70, ..., 2.00]

H_batch.shape = torch.Size([16, 200, 200])

[4.9342, 4.9344, 4.9345, ..., 4.9347]  ← all ≈ pi^2 / 2

1import torch
2import math
3
4device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
5
6# -------- 1.  Setup ------------------------------------------------------
7L  = 1.0
8N  = 200
9dx = L / (N + 1)
10
11# Build the same tridiagonal Hamiltonian, but as a torch.Tensor
12diag = torch.full((N,),      1.0 / dx**2,       device=device)
13off  = torch.full((N - 1,), -0.5 / dx**2,       device=device)
14H = torch.diag(diag) + torch.diag(off, 1) + torch.diag(off, -1)
15
16# -------- 2.  Diagonalize ------------------------------------------------
17# torch.linalg.eigh is the GPU/autograd-friendly analogue of np.linalg.eigh.
18# It accepts batched matrices, so a stack of Hamiltonians can be solved
19# in one call -- very useful when sweeping over L or V parameters.
20energies, states = torch.linalg.eigh(H)
21
22# -------- 3.  Compare to the exact result --------------------------------
23print(f"{'n':>3}  {'numerical':>12}  {'exact':>12}  {'error':>10}")
24for n in range(1, 6):
25    exact = (n * math.pi) ** 2 / (2 * L**2)
26    num   = energies[n - 1].item()
27    print(f"{n:3d}  {num:12.4f}  {exact:12.4f}  {abs(num - exact):10.2e}")
28
29# -------- 4.  Sweep over many box lengths in one shot --------------------
30# Vary L and watch E_1 scale like 1 / L^2.
31Ls    = torch.linspace(0.5, 2.0, 16, device=device)
32batch = []
33for L_val in Ls:
34    dx = (L_val / (N + 1)).item()
35    d  = torch.full((N,),      1.0 / dx**2, device=device)
36    o  = torch.full((N - 1,), -0.5 / dx**2, device=device)
37    batch.append(torch.diag(d) + torch.diag(o, 1) + torch.diag(o, -1))
38H_batch = torch.stack(batch)            # shape (16, N, N)
39
40E_batch, _ = torch.linalg.eigh(H_batch) # shape (16, N)
41E1_curve = E_batch[:, 0]                # ground-state energy for each L
42
43# E_1 * L^2 should be (very close to) the constant pi^2 / 2 ~ 4.9348
44print("E_1 * L^2 across the sweep (should be constant ~= pi^2/2):")
45print((E1_curve * Ls**2).cpu().numpy().round(4))

Why This Toy Model Matters

1. Quantum dots — colour by size

Real semiconductor quantum dots are three-dimensional infinite-well cousins of the model we just solved. Their emission wavelength is controlled almost entirely by their physical size: $E \propto 1/L^2 \Rightarrow \lambda \propto L^2$. Grow a CdSe dot at 2 nm and it glows blue; grow it at 6 nm and it glows red. Same chemistry, different box.

2. Quantum wells in lasers

Inside a semiconductor laser is a thin layer of low-bandgap material sandwiched between higher-bandgap layers — a near-perfect one-dimensional box for charge carriers. The quantized energy levels we just derived determine the laser's emission frequency.

3. Conjugated molecules and the free-electron model

In molecules like β-carotene, π-electrons are roughly free along the carbon backbone. Treating that backbone as an infinite well predicts the absorption colour of the molecule — it is why carrots are orange and tomatoes are red.

4. Numerical PDE solvers and ML for physics

Discretizing a Hamiltonian and diagonalizing it is the simplest case of the technique behind density-functional theory, neural quantum states, and physics-informed neural networks. The matrix we just diagonalized is exactly the building block those methods scale up to.

Test Your Understanding

Summary

Key Takeaways

1. Quantization is geometric. The discreteness of $E_n$ comes from a boundary condition, not from any quantum-mechanical postulate beyond Schrödinger's equation itself.
2. Localizing costs energy. The ground state has strictly positive energy. There is no way to put a particle in a box of finite size with zero kinetic energy.
3. Eigenstates are stationary; only superpositions move. The clock rate of the motion is set by the energy gap.
4. The orthonormal basis ψ_n is a Fourier sine basis. Any state on $[0,L]$ vanishing at the walls expands into them, and time evolution becomes phase-rotation of the expansion coefficients.
5. Numerics confirm everything. A 200-point finite-difference Hamiltonian recovers the analytic spectrum to four decimals, in both NumPy and PyTorch, and a single batched torch.linalg.eigh sweep proves $E_1 \propto 1/L^2$ directly.

Coming Next: In the next section we replace the rectangular potential with a parabola $V(x) = \tfrac{1}{2}m\omega^2 x^2$ and meet the quantum harmonic oscillator — the second pillar of every quantum textbook. There the eigenfunctions are Hermite polynomials times a Gaussian, the energies become evenly spaced $E_n = \hbar\omega(n+\tfrac{1}{2})$, and we discover the famous $\tfrac{1}{2}\hbar\omega$ zero-point energy of empty space.