Simulation of a quantum particle in a gravitational field.

A full simulation of a quantum particle in a gravitational field. 

In this project, we do three things:

1. Stationary State Analysis:

We set up the time‐independent Schrödinger equation for a particle in a gravitational potential

-\frac{\hbar^2}{2m}\frac{d^2\psi}{dz^2} + mg\,z\,\psi = E\,\psi,

2. Analytical Comparison via Airy Functions:

For the linear potential with a hard wall at , the ground state solution is known (up to normalization) to be given by an Airy function

\psi(z) \propto \mathrm{Ai}\Bigl(\frac{z-E/(mg)}{z_0}\Bigr),

z_0=\Bigl(\frac{\hbar^2}{2m^2g}\Bigr)^{1/3}.

3. Time Evolution Simulation:

We then propagate an initial Gaussian wave packet using the Crank–Nicolson scheme to solve

i\hbar\frac{\partial\psi}{\partial t} = H\psi,

---

%% Quantum Particle in a Gravitational Field: Full Simulation and Analysis

% This script performs:

% (i) A stationary state analysis of a quantum particle in a gravitational field,

% (ii) A comparison of the numerical ground state with the analytical Airy function solution,

% (iii) A time-dependent simulation of an initial Gaussian wave packet using the Crank–Nicolson method.

%

% The physical system obeys:

% - (ħ^2/(2m)) d^2ψ/dz^2 + mg*z*ψ = E ψ, with ψ(0)=ψ(L)=0.

%

% The analytical ground state is approximated by:

% ψ(z) ∝ Ai((z - E/(mg)) / z0),

% where z0 = (ħ^2/(2*m^2*g))^(1/3) and the first zero of Ai(–x) is ~2.3381.

%

% This code is entirely original while drawing on standard methods from quantum mechanics.

clear; clc; close all;

%% 1. Define Physical Constants and Domain

hbar = 1.0545718e-34; % Reduced Planck constant [J·s]

m = 1.675e-27; % Mass of neutron [kg]

g = 9.81; % Gravitational acceleration [m/s^2]

% Characteristic length scale for gravitational quantum states

z0 = (hbar^2/(2*m^2*g))^(1/3);

% Spatial domain

L = 100e-6; % Domain length: 100 micrometers

N = 1500; % Number of grid points

z = linspace(0, L, N)'; % Spatial grid (column vector)

dz = z(2) - z(1); % Spatial step size

%% 2. Build the Hamiltonian via Finite Differences

% Second derivative using central finite differences:

D2 = (diag(ones(N-1,1), 1) - 2*diag(ones(N,1), 0) + diag(ones(N-1,1), -1)) / dz^2;

% Gravitational potential: V(z) = m*g*z

V = m * g * z;

H = - (hbar^2/(2*m)) * D2 + diag(V);

%% 3. Stationary States: Eigenvalue Problem

numStates = 4; % Number of eigenstates to compute

[eigVecs, eigValMat] = eigs(H, numStates, 'smallestreal');

eigValues = diag(eigValMat);

[sortedEigValues, sortIdx] = sort(eigValues);

stationaryStates = eigVecs(:, sortIdx);

% Normalize the eigenfunctions

for k = 1:numStates

    stationaryStates(:,k) = stationaryStates(:,k) / sqrt(trapz(z, abs(stationaryStates(:,k)).^2));

end

% Display the computed eigenvalues in electron-volts (eV)

energies_eV = sortedEigValues / 1.60218e-19;

disp('Computed Energy Levels (eV):');

disp(energies_eV);

% Plot the probability densities of the stationary states

figure;

hold on;

colors = jet(numStates);

for k = 1:numStates

    plot(z*1e6, abs(stationaryStates(:,k)).^2, 'Color', colors(k,:), 'LineWidth', 2);

end

xlabel('Position (micrometers)');

ylabel('Probability Density');

title('Stationary States in a Gravitational Field');

legend('State 1','State 2','State 3','State 4');

grid on;

hold off;

%% 4. Analytical Ground State via Airy Function

% For a linear potential with a hard wall at z=0, the ground state energy is determined by:

% Ai(-E/(m*g*z0)) = 0.

% The first zero of the Airy function is approximately a1 = 2.3381, hence:

a1 = 2.3381;

E1_analytical = m * g * a1 * z0; % Ground state energy

% The analytical ground state (unnormalized) is:

psi_air = airy((z - E1_analytical/(m*g)) / z0);

% Normalize the analytical solution:

psi_air = psi_air / sqrt(trapz(z, psi_air.^2));

% Plot the numerical ground state vs. the analytical Airy function solution:

figure;

plot(z*1e6, abs(stationaryStates(:,1)).^2, 'b-', 'LineWidth', 2); hold on;

plot(z*1e6, psi_air.^2, 'r--', 'LineWidth', 2);

xlabel('Position (micrometers)');

ylabel('Probability Density');

title('Ground State: Numerical vs. Analytical Airy Function');

legend('Numerical Ground State','Analytical Airy Function');

grid on;

%% 5. Time Evolution using Crank–Nicolson

% Define an initial Gaussian wave packet:

z_center = 50e-6; % Center of the packet at 50 micrometers

sigma = 5e-6; % Packet width (standard deviation)

psi_initial = exp(-((z - z_center).^2) / (2*sigma^2));

psi_initial = psi_initial / sqrt(trapz(z, abs(psi_initial).^2)); % Normalize

% Time evolution parameters:

dt = 1e-7; % Time step [s]

totalTime = 1e-4; % Total simulation time [s]

numTimeSteps = round(totalTime / dt);

% Precompute Crank–Nicolson matrices:

I = eye(N);

A = I - 1i * dt/(2*hbar) * H;

B = I + 1i * dt/(2*hbar) * H;

% Initialize wavefunction for time evolution:

psi_t = psi_initial;

% Set up figure for animation:

figure;

hPlot = plot(z*1e6, abs(psi_t).^2, 'LineWidth', 2);

xlabel('Position (micrometers)');

ylabel('Probability Density');

title('Time Evolution of a Quantum Wave Packet');

axis([0 L*1e6 0 max(abs(psi_t).^2)*1.2]);

grid on;

% Time-stepping loop:

for t = 1:numTimeSteps

    psi_t = A \ (B * psi_t); % Crank–Nicolson update

    psi_t = psi_t / sqrt(trapz(z, abs(psi_t).^2)); % Renormalize

    

    % Update the plot every 20 time steps for performance:

    if mod(t,20) == 0

        set(hPlot, 'YData', abs(psi_t).^2);

        title(sprintf('Time Evolution: t = %.2e s', t*dt));

        drawnow;

    end

end

---

Detailed Explanation

1. Parameter Setup:

The script begins by defining fundamental constants (ħ, m, g) and calculates the characteristic length scale which sets the scale for gravitational quantum states. The spatial domain is defined from 0 to 100 micrometers, discretized into 1500 points.

2. Hamiltonian Construction:

The second derivative is approximated using a central finite-difference scheme. The Hamiltonian operator is built as

H = -\frac{\hbar^2}{2m} \frac{d^2}{dz^2} + mg\,z,

3. Stationary States:

We solve the eigenvalue problem using MATLAB’s sparse eigenvalue solver. The lowest four eigenstates are computed and normalized. Their corresponding energies are displayed (converted into electron-volts) and plotted.

4. Analytical Comparison:

Using the known result for a particle in a linear potential, the analytical ground state is expressed in terms of the Airy function. With the first zero , the ground state energy is given by

E_1 = m\,g\,(a_1\,z_0).

5. Time Evolution:

An initial Gaussian wave packet is defined and propagated in time using the Crank–Nicolson method. This implicit method is unconditionally stable and updates the wave function via a linear system at each time step. The evolving probability density is animated to illustrate the dynamics under the gravitational potential.