Quantum Mechanics Using MATLAB
Quantum mechanics is a fundamental branch of physics that deals with the behavior of matter and light on the atomic and subatomic scales. MATLAB provides robust tools and functions to model and solve quantum mechanical problems, such as solving the Schrödinger equation, simulating quantum systems, and analyzing eigenvalues and eigenvectors.
Solving Baiscs of Quantum Mechanics with MATLAB
1. Representing Quantum States
Quantum states are represented as vectors or matrices in MATLAB. For instance, the state vector of a two-level quantum system (qubit) can be defined as:
% Representing a quantum state
psi = [1; 0]; % Column vector representing the |0> state
disp(psi);
1
0
Similarly, the |1> state is represented as:
psi_1 = [0; 1]; % Column vector representing the |1> state
disp(psi_1);
0
1
2. Solving Schrödinger Equation
The time-independent Schrödinger equation is given by:
Hψ = Eψ, where H is the Hamiltonian operator, ψ is the wave function, and E is the energy eigenvalue.
\[ H \psi = E \psi \]
To solve this equation in MATLAB:
% Define the Hamiltonian matrix
H = [2 -1; -1 2];
% Solve for eigenvalues and eigenvectors
[E, psi] = eig(H);
disp('Eigenvalues:');
disp(diag(E));
disp('Eigenvectors:');
disp(psi);
Eigenvalues:
1.0000
3.0000
Eigenvectors:
-0.7071 0.7071
0.7071 0.7071
3. Particle in a Box
The particle-in-a-box model is one of the simplest quantum mechanical systems. The eigenfunctions are:
ψn(x) = sqrt(2/L) * sin(nπx/L), where L is the box length.
\[ \psi_n(x) = \sqrt{\frac{2}{L}} \sin\left(\frac{n\pi x}{L}\right) \]
To calculate and plot these eigenfunctions in MATLAB:
% Parameters
L = 1; % Box length
n = 1:3; % First three states
x = linspace(0, L, 100); % Position array
% Calculate and plot eigenfunctions
for k = n
psi_n = sqrt(2/L) * sin(k * pi * x / L);
plot(x, psi_n); hold on;
end
title('Particle in a Box Eigenfunctions');
xlabel('x'); ylabel('\psi_n(x)');
legend('n=1', 'n=2', 'n=3');
hold off;
4. Quantum Harmonic Oscillator
The quantum harmonic oscillator is described by the Schrödinger equation:
Hψ = (p²/2m + 1/2 kx²)ψ = Eψ
\[ H \psi = \left( \frac{p^2}{2m} + \frac{1}{2} kx^2 \right) \psi = E \psi \]
Here is a MATLAB code to numerically solve the oscillator:
% Discretize space
x = linspace(-5, 5, 200);
dx = x(2) - x(1);
% Kinetic and potential energy
T = -1/(2*dx^2) * diag(ones(199,1), 1) + diag(-2*ones(200,1)) + diag(ones(199,1), -1);
V = diag(0.5 * x.^2);
% Hamiltonian
H = T + V;
% Eigenvalues and eigenfunctions
[E, psi] = eig(H);
disp('Ground state energy:');
disp(E(1,1));
Advanced Quantum Mechanics Problems Using MATLAB
Example 1: Time Evolution of a Wave Packet
The time evolution of a wave packet is governed by the time-dependent Schrödinger equation. In this example, we simulate the time evolution of a Gaussian wave packet in a free particle system:
% Parameters
L = 10; % Length of the domain
N = 200; % Number of points
x = linspace(-L/2, L/2, N); % Spatial domain
dx = x(2) - x(1);
k0 = 5; % Initial wave vector
sigma = 0.5; % Width of the wave packet
t = 0:0.01:2; % Time array
% Initial Gaussian wave packet
psi0 = exp(-x.^2/(2*sigma^2)) .* exp(1i*k0*x);
psi0 = psi0 / sqrt(trapz(x, abs(psi0).^2)); % Normalize
% Kinetic energy operator in Fourier space
k = linspace(-pi/dx, pi/dx, N);
K = fftshift(k.^2 / 2);
% Time evolution
psi_t = zeros(N, length(t));
psi_t(:, 1) = psi0;
for j = 2:length(t)
psi_k = fft(psi_t(:, j-1));
psi_k = psi_k .* exp(-1i * K * (t(j) - t(j-1)));
psi_t(:, j) = ifft(psi_k);
end
% Visualization
for j = 1:length(t)
plot(x, abs(psi_t(:, j)).^2);
title(['Time Evolution at t = ', num2str(t(j))]);
xlabel('x'); ylabel('|\psi(x, t)|^2');
drawnow;
end
Example 2: Quantum Tunneling Through a Barrier
In quantum mechanics, particles can tunnel through potential barriers even if their energy is less than the barrier height. This example demonstrates tunneling through a Gaussian barrier:
% Parameters
L = 10; % Domain length
N = 200; % Number of grid points
x = linspace(-L/2, L/2, N);
dx = x(2) - x(1);
V0 = 5; % Barrier height
sigma = 0.5; % Barrier width
k0 = 3; % Wave vector
% Potential barrier
V = V0 * exp(-x.^2 / (2 * sigma^2));
% Hamiltonian matrix
T = -1/(2*dx^2) * (diag(ones(N-1, 1), 1) + diag(ones(N-1, 1), -1) - 2*diag(ones(N, 1)));
H = T + diag(V);
% Solve eigenvalues and eigenfunctions
[E, psi] = eig(H);
disp('Eigenvalues:');
disp(diag(E));
% Plot potential and eigenfunctions
plot(x, V, 'k', 'LineWidth', 2); hold on;
for n = 1:5
plot(x, abs(psi(:, n)).^2 + diag(E(n, n)));
end
xlabel('x'); ylabel('Potential and Eigenfunctions');
hold off;
Example 3: Variational Method for Ground State Energy
The variational method is a powerful technique to approximate the ground state energy of a quantum system. Here, we use a trial wavefunction to estimate the ground state energy of a quantum harmonic oscillator:
% Trial wavefunction parameters
alpha = 1; % Variational parameter
L = 5; % Domain length
N = 200; % Number of grid points
x = linspace(-L/2, L/2, N);
dx = x(2) - x(1);
% Trial wavefunction
psi_trial = exp(-alpha * x.^2);
psi_trial = psi_trial / sqrt(trapz(x, abs(psi_trial).^2)); % Normalize
% Kinetic and potential energy
T = -0.5 * diff([0; psi_trial]) ./ dx;
T = T(1:end-1) .* psi_trial(1:end-1);
V = 0.5 * x.^2 .* abs(psi_trial).^2;
% Total energy
E_trial = trapz(x(1:end-1), T + V);
disp('Estimated Ground State Energy:');
disp(E_trial);
Example 4: Quantum Entanglement and Bell’s Inequality
Quantum entanglement can be tested using Bell’s inequality. Here, we simulate the CHSH inequality:
% Bell's inequality simulation
theta_A = 0; theta_A_prime = pi/4;
theta_B = pi/8; theta_B_prime = 3*pi/8;
% Measurement correlation function
E = @(theta1, theta2) cos(2 * (theta1 - theta2));
% Calculate CHSH inequality
S = abs(E(theta_A, theta_B) - E(theta_A, theta_B_prime) + ...
E(theta_A_prime, theta_B) + E(theta_A_prime, theta_B_prime));
disp('CHSH Inequality Value (S):');
disp(S);
% Verify violation
if S > 2
disp('Violation of Bell''s inequality detected.');
else
disp('No violation detected.');
end
Example 5: Solving the Schrödinger Equation for a Quantum Harmonic Oscillator
The Schrödinger equation for a quantum harmonic oscillator has solutions in terms of Hermite polynomials. Here, we numerically solve the equation using finite difference methods:
% Parameters
N = 1000; % Number of grid points
L = 10; % Length of the domain
x = linspace(-L/2, L/2, N); % Grid points
dx = x(2) - x(1); % Step size
m = 1; % Mass of the particle
hbar = 1; % Reduced Planck's constant
omega = 1; % Angular frequency
% Potential energy: Harmonic potential
V = 0.5 * m * omega^2 * x.^2;
% Hamiltonian matrix using finite difference
T = -hbar^2 / (2 * m * dx^2); % Kinetic energy coefficient
H = diag(2 * T + V) + diag(-T * ones(N-1, 1), 1) + diag(-T * ones(N-1, 1), -1);
% Eigenvalues and eigenfunctions
[psi, E] = eig(H);
E = diag(E); % Extract eigenvalues
% Plot the first few wavefunctions
figure;
hold on;
plot(x, psi(:, 1).^2, 'DisplayName', 'Ground State');
plot(x, psi(:, 2).^2, 'DisplayName', 'First Excited State');
plot(x, psi(:, 3).^2, 'DisplayName', 'Second Excited State');
legend;
xlabel('x');
ylabel('|\psi(x)|^2');
title('Wavefunctions of the Quantum Harmonic Oscillator');
Example 6: Time Evolution of a Gaussian Wave Packet
We study the time evolution of a Gaussian wave packet in free space by numerically solving the time-dependent Schrödinger equation.
% Parameters
N = 500; % Number of grid points
L = 20; % Length of the domain
x = linspace(-L/2, L/2, N);
dx = x(2) - x(1);
hbar = 1; % Reduced Planck's constant
m = 1; % Mass
% Initial wave packet
x0 = -5; % Initial position
k0 = 5; % Initial wave vector
sigma = 1; % Width of the wave packet
psi0 = exp(-(x - x0).^2 / (2 * sigma^2)) .* exp(1i * k0 * x);
psi0 = psi0 / sqrt(sum(abs(psi0).^2) * dx); % Normalize
% Time evolution operator
T = -hbar^2 / (2 * m * dx^2); % Kinetic energy coefficient
H = diag(-2 * T * ones(N, 1)) + diag(T * ones(N-1, 1), 1) + diag(T * ones(N-1, 1), -1);
U = expm(-1i * H * 0.01 / hbar); % Time evolution matrix for small dt
% Time evolution
psi = psi0;
figure;
for t = 1:200
psi = U * psi;
if mod(t, 20) == 0 % Plot every 20 steps
plot(x, abs(psi).^2, 'DisplayName', ['t = ' num2str(t * 0.01)]);
hold on;
end
end
xlabel('x');
ylabel('|\psi(x)|^2');
title('Time Evolution of a Gaussian Wave Packet');
legend;
Example 7: Solving the Hydrogen Atom Radial Equation
We numerically solve the radial Schrödinger equation for the hydrogen atom, obtaining energy levels and wavefunctions.
% Parameters
N = 1000; % Number of grid points
r_max = 20; % Maximum radial distance
r = linspace(1e-5, r_max, N); % Radial grid (avoid r=0 for singularity)
dr = r(2) - r(1);
Z = 1; % Atomic number (hydrogen)
m = 1; % Reduced mass
hbar = 1; % Reduced Planck's constant
% Effective potential
V = -Z ./ r + hbar^2 ./ (2 * m * r.^2);
% Hamiltonian matrix
T = -hbar^2 / (2 * m * dr^2);
H = diag(2 * T + V) + diag(-T * ones(N-1, 1), 1) + diag(-T * ones(N-1, 1), -1);
% Eigenvalues and eigenfunctions
[psi, E] = eig(H);
E = diag(E); % Extract eigenvalues
% Plot the radial wavefunctions
figure;
plot(r, psi(:, 1).^2, 'DisplayName', 'Ground State');
hold on;
plot(r, psi(:, 2).^2, 'DisplayName', 'First Excited State');
plot(r, psi(:, 3).^2, 'DisplayName', 'Second Excited State');
legend;
xlabel('r');
ylabel('|\psi(r)|^2');
title('Radial Wavefunctions of Hydrogen Atom');
Useful MATLAB Functions for Quantum Mechanics
Practice Questions
Test Yourself
1. Solve the Schrödinger equation for a particle in a finite square well potential.
2. Simulate the time evolution of a wave packet in a potential well.
3. Compute and visualize the eigenfunctions of a quantum harmonic oscillator.